bhm.mcmc {mederrRank}

Markov Chain Monte Carlo Estimation (Step 1) of the Bayesian Hierarchical Model for Identifying the Most Harmful Medication Errors

Description

This function implements the Markov Chain Monte Carlo estimation methodology for the Bayesian hierarchical model described in Myers et al. (2011).

Usage

bhm.mcmc(dat, nsim = 2000, burnin = 500, scale.factor = 1,
	adaptive.int = 100, adaptive.max = 1000, prior = NULL,
	init = NULL, tuneD = NULL, tuneT = NULL)

Arguments

dat

an object of class "mederrData".

nsim

number of iterations.

burnin

number of burn-in iterations.

scale.factor

scale factor of the random effects proposal distribution.

adaptive.int

iteration interval at which the standard error of the random effects proposal distribution is updated.

adaptive.max

last iteration at which the standard error of the random effects proposal distribution is updated.

prior

an optional list of the hyperparameters values; see the Details section below.

init

an optional list of initial values for the model parameters; see the Details section below.

tuneD

an optional vector of the \delta_j proposal distribution variances.

tuneT

an optional vector of the \theta_i proposal distribution variances.

Details

The Bayesian hierarchical model (with crossed random effects) implemented here for identifying the medication error profiles with the largest log odds of harm is

y_{ij} | N_{ij}, p_{ij} \sim Bin(N_{ij},p_{ij})

logit(p_{ij}) = \gamma + \theta_i + \delta_j

\theta_i | \sigma, \eta, k \sim St(0,\sigma,k,\eta), \qquad i=1,\ldots,n

\delta_j | \tau^2 \sim N(0,\tau^2), \qquad j=1,\ldots,J

\gamma \sim N(g,G)

\sigma^2 \sim IG(a_1,b_1)

\tau^2 \sim IG(a_2,b_2)

k \sim Unif(0,\infty)

\eta \sim Unif(0,\infty),

where N_{ij} denotes the number of times that the error profile i is cited on a report from hospital j and y_{ij} is the corresponding number of times that profile i in hospital j was reported with harm. This function implements the first model estimation step in which the values k = \infty and k = 1, i.e. a symmetric normal distribution, is forced for the error profiles' random effects. A sample from the joint posterior distribution of all other parameters via Markov Chain Monte Carlo with adaptive Metropolis steps for each set of random effects is obtained. For more details see Myers et al. (2011).

Value

bhm.mcmc returns an object of the class "mederrFit".

Author(s)

Sergio Venturini [e-mail],

Jessica A. Myers [e-mail]

References

Myers, J. A., Venturini, S., Dominici, F. and Morlock, L. (2011), "Random Effects Models for Identifying the Most Harmful Medication Errors in a Large, Voluntary Reporting Database". Technical Report.

See Also

bhm.resample, mederrData, mederrFit.

Examples

## Not run: 
data("simdata", package = "mederrRank")
summary(simdata)

fit <- bhm.mcmc(simdata, nsim = 1000, burnin = 500, scale.factor = 1.1)
resamp <- bhm.resample(fit, simdata, p.resample = .1,
	k = c(3, 6, 10, 30, 60, Inf), eta = c(.5, .8, 1, 1.25, 2))
fit2 <- bhm.constr.resamp(fit, resamp, k = 3, eta = .8)
plot(fit, fit2, simdata)

theta0 <- c(10, 6, 100, 100, .1)
ans <- mixnegbinom.em(simdata, theta0, 50000, 0.01, se = TRUE,
	stratified = TRUE)

summary(fit2, ans, simdata)
## End(Not run)