BMD {CoMiRe}

Benchmark dose

Description

Benchmark dose associated to a particular risk

Usage

BMD(level, risk, x, alpha=0.05)

Arguments

level

dose level of interest.

risk

summary.risk$mcmc.risk from the output of add.risk function.

x

numeric vector for the covariate relative to the dose of exposure used in comire.gibbs.

alpha

level of the credible bands.

Value

A dataframe containing as variables:

• q the dose level of interest.

• BMD the benchmark dose.

• low lower credible limit.

• upp upper credible limit.

• BMDL a more conservative benchmark dose.

Author(s)

Antonio Canale

Examples

{
data(CPP)
attach(CPP)

n <- NROW(CPP)
J <- H <- 10

premature <- as.numeric(gestage<=37)

mcmc <- list(nrep=5000, nb=2000, thin=5, ndisplay=4)

## too few iterations to be meaningful. see below for safer and more comprehensive results

mcmc <- list(nrep=10, nb=2, thin=1, ndisplay=4) 

prior <- list(mu.theta=mean(gestage), k.theta=10, eta=rep(1, J)/J, 
              alpha=rep(1,H)/H, a=2, b=2, J=J, H=H)
              
fit.dummy <- comire.gibbs(gestage, dde, family="continuous", 
                     mcmc=mcmc, prior=prior, seed=1, max.x=180)
                     
risk.data <- add.risk(y = gestage, x = dde, fit = fit.dummy, mcmc = mcmc, 
    a = 37, x.grid = seq(0, max(dde), length = 100))
bmd.data <- BMD(seq(0,.20, length=50), risk.data$mcmc.risk, 
x=seq(0,max(dde), length=100), alpha=0.05)
bmd.plot(bmd.data)       
                     

## safer procedure with more iterations (it may take some time)

mcmc <- list(nrep=5000, nb=2000, thin=5, ndisplay=4)

## Fit the model for continuous y 

prior <- list(mu.theta=mean(gestage), k.theta=10, eta=rep(1, J)/J, 
              alpha=rep(1,H)/H, a=2, b=2, J=J, H=H)
              
fit <- comire.gibbs(gestage, dde, family="continuous", 
                     mcmc=mcmc, prior=prior, seed=5, max.x=180)


risk.data <- add.risk(y = gestage, x = dde, fit = fit, mcmc = mcmc,
a = 37, x.grid = seq(0, max(dde), length = 100))
bmd.data <- BMD(seq(0,.20, length=50), risk.data$mcmc.risk, 
x=seq(0,max(dde), length=100), alpha=0.05)
bmd.plot(bmd.data)       


}