{
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))
riskplot(risk.data$summary.risk, xlab="DDE", x = dde, xlim = c(0,150))
# \donttest{
## 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))
riskplot(risk.data$summary.risk, xlab="DDE",
x = dde, xlim = c(0,150))
# }
}
Run the code above in your browser using DataLab