Title: | Bayesian Methods for Identifying the Most Harmful Medication Errors |
---|---|
Description: | Two distinct but related statistical approaches to the problem of identifying the combinations of medication error characteristics that are more likely to result in harm are implemented in this package: 1) a Bayesian hierarchical model with optimal Bayesian ranking on the log odds of harm, and 2) an empirical Bayes model that estimates the ratio of the observed count of harm to the count that would be expected if error characteristics and harm were independent. In addition, for the Bayesian hierarchical model, the package provides functions to assess the sensitivity of results to different specifications of the random effects distributions. |
Authors: | Sergio Venturini, Jessica Myers |
Maintainer: | Sergio Venturini <[email protected]> |
License: | GPL (>= 2) | file LICENSE |
Version: | 0.1.0 |
Built: | 2024-12-05 06:53:47 UTC |
Source: | CRAN |
Two distinct but related statistical approaches to the problem of identifying the combinations of medication error characteristics that are more likely to result in harm are implemented in this package: 1) a Bayesian hierarchical model with optimal Bayesian ranking on the log odds of harm, and 2) an empirical Bayes model that estimates the ratio of the observed count of harm to the count that would be expected if error characteristics and harm were independent. In addition, for the Bayesian hierarchical model, the package provides functions to assess the sensitivity of results to different specifications of the random effects distributions.
The package is loaded with the usual library(mederrRank)
command. The most important functions are bhm.mcmc
, bhm.resample
and mixnegbinom.em
.
Sergio Venturini, Jessica Myers
Maintainer: Sergio Venturini <[email protected]>
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.
bayes.rank
,
bhm.mcmc
,
bhm.resample
,
mixnegbinom.em
.
This function estimates the ranks of the log odds of harm of the various medication error profiles as described in Myers et al. (2011).
bayes.rank(model)
bayes.rank(model)
model |
a mederrFit object. |
Using the posterior samples of the , the function estimates the ranks of the log odds of harm of the various error profiles. Optimal Bayesian ranking gives estimates of rank for profile
as
where is the posterior probability that
.
bayes.rank
returns the numerical vector of Optimal Bayesian ranks for the chosen mederrFit model (see the references for the details).
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
## Not run: data("simdata", package = "mederrRank") summary(simdata) fit <- bhm.mcmc(simdata, nsim = 1000, burnin = 500, scale.factor = 1.1) ranks <- bayes.rank(fit) summary(ranks) ## End(Not run)
## Not run: data("simdata", package = "mederrRank") summary(simdata) fit <- bhm.mcmc(simdata, nsim = 1000, burnin = 500, scale.factor = 1.1) ranks <- bayes.rank(fit) summary(ranks) ## End(Not run)
This function represents the "constructor" function for the resampling procedure used in this package. bhm.resample
calculates the importance ratios, and performs the sampling, and then this function constructs the resampled model based on that information.
bhm.constr.resamp(model, resample, k, eta)
bhm.constr.resamp(model, resample, k, eta)
model |
|
resample |
an object of class "mederrResample". |
k |
|
eta |
|
Deviations from the normal, i.e. , random effects distribution using a different pair of
and
values are considered. The methodology implemented here is the importance link function resampling approach introduce by MacEachern and Peruggia (2000): based on the
chain, new posterior samples under a new set of
values is obtained.
bhm.constr.resamp
returns an object of the class "mederrFit".
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
MacEachern, S. and Peruggia, M. (2000), "Importance Link Function Estimation for Markov Chain Monte Carlo Methods", Journal of Computational and Graphical Statistics, 9, 99-121.
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.
bhm.mcmc
,
bhm.resample
,
mederrData
,
mederrFit
.
## 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)
## 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)
This function implements the Markov Chain Monte Carlo estimation methodology for the Bayesian hierarchical model described in Myers et al. (2011).
bhm.mcmc(dat, nsim = 2000, burnin = 500, scale.factor = 1, adaptive.int = 100, adaptive.max = 1000, prior = NULL, init = NULL, tuneD = NULL, tuneT = NULL)
bhm.mcmc(dat, nsim = 2000, burnin = 500, scale.factor = 1, adaptive.int = 100, adaptive.max = 1000, prior = NULL, init = NULL, tuneD = NULL, tuneT = NULL)
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 |
tuneT |
an optional vector of the |
The Bayesian hierarchical model (with crossed random effects) implemented here for identifying the medication error profiles with the largest log odds of harm is
where denotes the number of times that the error profile
is cited on a report from hospital j and
is the corresponding number of times that profile
in hospital
was reported with harm.
This function implements the first model estimation step in which the values
and
, 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).
bhm.mcmc
returns an object of the class "mederrFit".
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.resample
,
mederrData
,
mederrFit
.
## 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)
## 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)
This function implements the transformation needed to apply the importance link function resampling methodology based on the Markov Chain Monte Carlo simulations obtained with the bhm.mcmc
command (see the References).
bhm.resample(model, dat, p.resample = 0.1, k, eta)
bhm.resample(model, dat, p.resample = 0.1, k, eta)
model |
|
dat |
an object of class "mederrData". |
p.resample |
proportion of simulations resampled from the |
k |
required vector of |
eta |
required vector of |
bhm.resample
returns an object of the class "mederrResample".
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
MacEachern, S. and Peruggia, M. (2000), "Importance Link Function Estimation for Markov Chain Monte Carlo Methods", Journal of Computational and Graphical Statistics, 9, 99-121.
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.
mederrData
,
mederrFit
,
bhm.mcmc
.
## 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)
## 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)
Density function for a mixture of two negative binomial distributions.
dmixnegbinom(x, theta, E, log.p = FALSE)
dmixnegbinom(x, theta, E, log.p = FALSE)
x |
vector of (non-negative integer) quantiles. |
theta |
vector of parameters for the negative binomial distribution mixture. |
E |
vector of (non-negative integer) expected counts. |
log.p |
logical; if TRUE, probabilities p are given as log(p). |
The mixture of two negative binomial distributions has density
where
for ,
and
.
The mixture of two negative binomial distributions represents the marginal distribution of the counts
coming from Poisson data with parameter
and a mixture of two gamma distributions as its prior. For details see the paper by Dumouchel (1999).
dmixnegbinom
gives the density corresponding to the E
and theta
values provided.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
## Not run: data("simdata", package = "mederrRank") ni <- simdata@numi theta0 <- c(10, 6, 100, 100, .1) ans <- mixnegbinom.em(simdata, theta0, 50000, 0.01, se = FALSE, stratified = TRUE) theta <- ans$theta.hat N.E <- cbind(ans$N[1:ni], ans$E[1:ni])[sort(ans$N[1:ni], index.return = TRUE)$ix, ] N.ix <- match(unique(N.E[, 1]), N.E[, 1]) N <- N.E[N.ix, 1] E <- N.E[N.ix, 2] dens <- dmixnegbinom(N, theta, E) hist(N.E[, 1], breaks = 40, freq = FALSE) points(N, dens) ## End(Not run)
## Not run: data("simdata", package = "mederrRank") ni <- simdata@numi theta0 <- c(10, 6, 100, 100, .1) ans <- mixnegbinom.em(simdata, theta0, 50000, 0.01, se = FALSE, stratified = TRUE) theta <- ans$theta.hat N.E <- cbind(ans$N[1:ni], ans$E[1:ni])[sort(ans$N[1:ni], index.return = TRUE)$ix, ] N.ix <- match(unique(N.E[, 1]), N.E[, 1]) N <- N.E[N.ix, 1] E <- N.E[N.ix, 2] dens <- dmixnegbinom(N, theta, E) hist(N.E[, 1], breaks = 40, freq = FALSE) points(N, dens) ## End(Not run)
Density function for the negative binomial distribution with parameters alpha
and prob
.
dnegbinom(x, alpha, prob, log.p = FALSE)
dnegbinom(x, alpha, prob, log.p = FALSE)
x |
vector of (non-negative integer) quantiles. |
alpha |
target for number of successful trials. Must be strictly positive, need not be integer. |
prob |
probability of success in each trial. |
log.p |
logical; if TRUE, probabilities p are given as log(p). |
The negative binomial distribution with parameters alpha
= and
prob
= has density
for and
.
This represents the number of failures which occur in a sequence of Bernoulli trials before a target number of successes is reached.
dnegbinom
gives the density corresponding to the alpha
and prob
values provided.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
Density function for the skewed t distribution with k
degrees of freedom, scale parameter sigma
and skewness eta
.
dst(x, sigma, k, eta)
dst(x, sigma, k, eta)
x |
vector of quantiles. |
sigma |
scale parameter ( |
k |
degrees of freedom ( |
eta |
skewness parameter ( |
This distribution is based on introducing skewing into the symmetric scaled t distribution, as described in Fernandez and Steel (1998).
The parameters characterizing the center (here set at 0) and the spread (sigma
) refer to the mean and standard deviation of the underlying symmetric distribution.
In the skewed t distribution, the centrality parameter defines the mode of the distribution, but it is no longer either the mean or the median. Similarly, in the skewed t distribution, sigma
still characterizes the spread, but it can no longer be interpreted directly as the standard deviation of the distribution.
dst
gives the density corresponding to the simga
, k
and eta
values provided.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
Fernandez, C. and Steel, M. (1998), "On Bayesian Modeling of Fat Tails and Skewness". Journal of the American Statistical Association, 93, 359-371.
Lee, K. and Thompson, S. (2008), "Flexible Parametric Models for Random-Effects Distributions". Statistics in Medicine, 27, 418-434.
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.
dt
.
This function computes the geometric mean of the empirical Bayes posterior distribution for the observed vs. expected count relative risk.
EBGM(eb.result)
EBGM(eb.result)
eb.result |
output of the |
For further details see DuMouchel (1999).
EBGM
returns the vector of geometric means.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
## 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) theta0 <- c(10, 6, 100, 100, .1) ans <- mixnegbinom.em(simdata, theta0, 50000, 0.01, se = FALSE, stratified = TRUE) ni <- simdata@numi rank(EBGM(ans)[1:ni]) summary(fit2, ans, simdata) ## End(Not run)
## 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) theta0 <- c(10, 6, 100, 100, .1) ans <- mixnegbinom.em(simdata, theta0, 50000, 0.01, se = FALSE, stratified = TRUE) ni <- simdata@numi rank(EBGM(ans)[1:ni]) summary(fit2, ans, simdata) ## End(Not run)
Parameters
This function computes the log-likelihood difference for the candidate random effects. It is a helper function and not meant to be used on its own.
llDiffD(dat, deltaj, cand, thetai, gamma, tau2)
llDiffD(dat, deltaj, cand, thetai, gamma, tau2)
dat |
data frame containing the observed sample counts. |
deltaj |
vector of previous accepted values for the |
cand |
vector of candidate values for the |
thetai |
vector of previous accepted values for the |
gamma |
last sampled value for the |
tau2 |
last sampled value for the |
For further details see Myers et al. (2011).
llDiffD
returns the vector of log-likelihood differences.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.constr.resamp
,
bhm.mcmc
,
bhm.resample
.
Parameters
This function computes the log-likelihood difference for the candidate random effects. It is a helper function and not meant to be used on its own.
llDiffT(dat, thetai, cand, deltaj, gamma, sigma2)
llDiffT(dat, thetai, cand, deltaj, gamma, sigma2)
dat |
data frame containing the observed sample counts. |
thetai |
vector of previous accepted values for the |
cand |
vector of candidate values for the |
deltaj |
vector of previous accepted values for the |
gamma |
last sampled value for the |
sigma2 |
last sampled value for the |
For further details see Myers et al. (2011).
llDiffT
returns the vector of log-likelihood differences.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.constr.resamp
,
bhm.mcmc
,
bhm.resample
.
This function computes the negative log-posterior distribution of the Bayesian hierarchical model described in Myers et al (2011). It is a helper function and not meant to be used on its own.
logp(theta, deltaj, sigma2, i, k, eta, dat)
logp(theta, deltaj, sigma2, i, k, eta, dat)
theta |
value of the error profile random effect at which the log.posterior distribution is calculated. |
deltaj |
vector of hospital random effect values. |
sigma2 |
scale parameter ( |
i |
error profile index for which the calculate of the log.posterior distribution is needed. |
k |
degrees of freedom ( |
eta |
skewness parameter ( |
dat |
an object of class "mederrData". |
For further details see Myers et al. (2011).
logp
returns a vector of log-posterior values.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.constr.resamp
,
bhm.mcmc
,
bhm.resample
.
and
This functions computes the unnormalized marginal posterior distributions for the and
parameters as described in Myers et al (2011).
logunpost(resample)
logunpost(resample)
resample |
an object of class "mederrResample". |
logunpost
is used in the plot method
for a mederrResample object.
logunpost
returns an array with the posterior distribution values.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.constr.resamp
,
bhm.mcmc
,
bhm.resample
,
mederrResample
.
This class encapsulates the data specification for a Bayesian Hierarchical Model used to identify the most harmful medication errors as described in Myers et al. (2011).
Objects can be created by calls of the form new("mederrData", data)
, where the data
argument has to be a matrix or a data frame object that contains the following (numeric) information for each error profile/hospital combination:
the number of times () that profile
in hospital
was reported with harm;
the total number of times () that the error profile
is cited on a report from hospital
,
the error profile identification code,
the hospital identification code.
data
:Object of class "data.frame"
; data in the standard data.frame
form.
size
:Object of class "numeric"
; total number of observations in the data set.
numi
:Object of class "numeric"
; number of error profiles available in the data set.
numj
:Object of class "numeric"
; number of hospitals available in the data set.
signature(x = "mederrData", y = "missing")
: Provides a pictorial representation for a sample of error profiles reported by some hospitals.
signature(object = "mederrData")
: Summarizes information about an mederrData
object.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bayes.rank
,
bhm.mcmc
,
bhm.resample
,
mixnegbinom.em
.
ng <- 50 i <- rep(1:ng, ng) j <- rep(1:ng, each = ng) N <- rpois(ng^2, 3 + .05*i - .01*j) + 1 theta_i <- rgamma(ng, 5, 5) - 4/5 delta_j <- rnorm(ng, 0, .2) logit <- -3 + theta_i[i] + delta_j[j] y <- rbinom(ng^2, N, exp(logit)/(1 + exp(logit))) simdata <- new("mederrData", data = cbind(y, N, i, j))
ng <- 50 i <- rep(1:ng, ng) j <- rep(1:ng, each = ng) N <- rpois(ng^2, 3 + .05*i - .01*j) + 1 theta_i <- rgamma(ng, 5, 5) - 4/5 delta_j <- rnorm(ng, 0, .2) logit <- -3 + theta_i[i] + delta_j[j] y <- rbinom(ng^2, N, exp(logit)/(1 + exp(logit))) simdata <- new("mederrData", data = cbind(y, N, i, j))
This class encapsulates the simulated Monte Carlo chains for the Bayesian Hierarchical Model as described in Myers et al. (2011) forcing a symmetric normal distribution on the ,
.
Objects can be created by calls of the form new("mederrFit", thetai, deltaj, gamma, sigma2, tau2, p.acc.i, p.acc.j, tune.theta, tune.delta, k, eta)
, but most often as the result of a call to bhm.mcmc
or to bhm.constr.resamp
.
thetai
:Object of class "matrix"
; simulated chains for the ,
, error profiles random effects; see
bhm.mcmc
.
deltaj
:Object of class "matrix"
; simulated chains for the ,
, hospitals random effects; see
bhm.mcmc
.
gamma
:Object of class "numeric"
; simulated chain for the parameter; see
bhm.mcmc
.
sigma2
:Object of class "numeric"
; simulated chain for the parameter; see
bhm.mcmc
.
tau2
:Object of class "numeric"
; simulated chain for the parameter; see
bhm.mcmc
.
p.acc.i
:Object of class "numeric"
; acceptance rates for the error profiles random effects.
p.acc.j
:Object of class "numeric"
; acceptance rates for the hospitals random effects.
tune.theta
:Object of class "numeric"
; last updated values of the working variances for the Metropolis step.
tune.delta
:Object of class "numeric"
; last updated values of the working variances for the Metropolis step.
k
:Object of class "numeric"
; value used in the simulation.
eta
:Object of class "numeric"
; value used in the simulation.
signature(x = "mederrFit", y = "mederrFit")
: Provides a graphical representation of the estimates.
signature(object = "mederrFit")
: Summarizes the information regarding the estimates.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bayes.rank
,
bhm.constr.resamp
,
bhm.mcmc
.
This class encapsulates the information needed to resample the Monte Carlo chains for the Bayesian Hierarchical Model as described in Myers et al. (2011) using user defined values for and
.
Objects can be created by calls of the form new("mederrResample", log.ir, samp, A, t.new, t.old, grd)
, but most often as the result of a call to bhm.resample
.
log.ir
:Object of class "array"
; logarithm of the importance ratio for each pair of values.
samp
:Object of class "array"
; resampled MCMC simulation indexes.
A
:Object of class "array"
; transformation ratio for each pair of values.
t.new
:Object of class "array"
; posterior modes using
.
t.old
:Object of class "numeric"
; posterior modes using user defined
values.
grd
:Object of class "list"
; grid of required values.
signature(x = "mederrResample", y = "missing")
: : Provides a graphical representation of a mederrResample
object.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bayes.rank
.
bhm.constr.resamp
,
bhm.mcmc
.
Subset of the MEDMARX data included for illustrative purposes only in the mederrRank
package.
data(MEDMARX)
data(MEDMARX)
An object of class mederrData
.
The data contained in this object are reproduced by gentle permission of Quantros, Inc., 690 N. McCarthy Blvd., Suite 200, Milpitas, CA 95035.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.mcmc
,
mederrData
,
mederrFit
.
data("MEDMARX", package = "mederrRank") summary(MEDMARX) plot(MEDMARX, nbins.err = 20, nbins.hosp = 10)
data("MEDMARX", package = "mederrRank") summary(MEDMARX) plot(MEDMARX, nbins.err = 20, nbins.hosp = 10)
This function provides the empirical Bayes estimates for the parameters theta
of a mixture of two negative binomial distributions (see dmixnegbinom
) using an Expectation-Maximization algorithm.
mixnegbinom.em(dat, theta0, maxiter = 50000, toler = 0.01, se = TRUE, stratified = FALSE)
mixnegbinom.em(dat, theta0, maxiter = 50000, toler = 0.01, se = TRUE, stratified = FALSE)
dat |
an object of class "mederrData". |
theta0 |
initial values for the parameters to be optimized over. |
maxiter |
a positive integer specifying the maximum number of iterations to be performed before the program is terminated. |
toler |
a positive scalar giving the tolerance at which the change in the log-likelihood is considered close enough to zero to terminate the algorithm. |
se |
logical; if TRUE the standard errors of the estimates are also returned. |
stratified |
logical; if TRUE the analysis will be performed by stratifying on the hospitals. |
For further details see Myers et al. (2011).
mixnegbinom.em
returns a list with components:
theta.hat |
The best set of parameters found. |
final.err |
The last change in the log-likelihood; it has to be smaller than the |
final.ll |
The likelihood value corresponding to |
final.score |
The log-likelihood score value corresponding to |
num.iter |
The number of iterations performed to find the proposed solution. |
se |
Only if argument |
N |
The vector of observed error profiles counts. |
E |
The vector of expected error profiles counts. |
prior |
A character string giving the prior used; for this function is set to "mixgamma", i.e. a mixture of two gamma distributions as in DuMouchel (1999). |
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
dmixnegbinom
,
EBGM
,
negbinom.em
.
## 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) ans$theta ans$se summary(fit2, ans, simdata) ## End(Not run)
## 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) ans$theta ans$se summary(fit2, ans, simdata) ## End(Not run)
This function computes the log-likelihood function for the mixture of two negative binomial distributions as described in dmixnegbinom
.
mixnegbinom.loglik(theta, N, E)
mixnegbinom.loglik(theta, N, E)
theta |
vector of parameter values. |
N |
vector of observed error profiles counts. |
E |
vector of expected error profiles counts. |
For further details see Myers et al. (2011).
mixnegbinom.loglik
returns the log-likelihood value for the negative binomial mixture.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
dmixnegbinom
,
mixnegbinom.em
,
mixnegbinom.score
.
This function computes the log-likelihood score for the mixture of two negative binomial distributions as described in dmixnegbinom
.
mixnegbinom.score(theta, N, E)
mixnegbinom.score(theta, N, E)
theta |
vector of parameter values. |
N |
vector of observed error profiles counts. |
E |
vector of expected error profiles counts. |
For further details see Myers et al. (2011).
mixnegbinom.score
returns the vector of log-likelihood score values for the negative binomial mixture.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
dmixnegbinom
,
mixnegbinom.em
,
mixnegbinom.loglik
.
This function provides the empirical Bayes estimates for the parameters theta
of a negative binomial distribution (see dnegbinom
) using an Expectation-Maximization algorithm.
negbinom.em(dat, theta0, maxiter = 50000, toler = 0.01, se = TRUE, stratified = FALSE)
negbinom.em(dat, theta0, maxiter = 50000, toler = 0.01, se = TRUE, stratified = FALSE)
dat |
an object of class "mederrData". |
theta0 |
initial values for the parameters to be optimized over. |
maxiter |
a positive integer specifying the maximum number of iterations to be performed before the program is terminated. |
toler |
a positive scalar giving the tolerance at which the scaled gradient is considered close enough to zero to terminate the algorithm. |
se |
logical; if TRUE the standard errors of the estimates are also returned. |
stratified |
logical; if TRUE the analysis will be performed by stratifying on the hospitals. |
For further details see Myers et al. (2011).
negbinom.em
returns a list with components:
theta.hat |
The best set of parameters found. |
final.err |
The last change in the log-likelihood; it has to be smaller than the |
final.ll |
The likelihood value corresponding to |
final.score |
The log-likelihood score value corresponding to |
num.iter |
The number of iterations performed to find the proposed solution. |
se |
Only if argument |
N |
The vector of observed error profiles counts. |
E |
The vector of expected error profiles counts. |
prior |
A character string giving the prior used; for this function is set to "gamma", i.e. a gamma distribution. |
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
dnegbinom
,
EBGM
,
mixnegbinom.em
.
data("simdata", package = "mederrRank") summary(simdata) ## Not run: 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) ## End(Not run) theta0 <- runif(2, 0, 5) ans <- negbinom.em(simdata, theta0, 50000, 0.01, se = TRUE, stratified = TRUE) ans$theta ans$se ## Not run: summary(fit2, ans, simdata) ## End(Not run)
data("simdata", package = "mederrRank") summary(simdata) ## Not run: 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) ## End(Not run) theta0 <- runif(2, 0, 5) ans <- negbinom.em(simdata, theta0, 50000, 0.01, se = TRUE, stratified = TRUE) ans$theta ans$se ## Not run: summary(fit2, ans, simdata) ## End(Not run)
This function computes the log-likelihood function for the mixture of two negative binomial distribution as described in dmixnegbinom
.
negbinom.loglik(theta, N, E)
negbinom.loglik(theta, N, E)
theta |
vector of parameter values. |
N |
vector of observed error profiles counts. |
E |
vector of expected error profiles counts. |
For further details see Myers et al. (2011).
negbinom.loglik
returns the log-likelihood value for the negative binomial distribution.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
dnegbinom
,
negbinom.em
,
negbinom.score
.
This function computes the log-likelihood score for the negative binomial distribution as described in dmixnegbinom
.
negbinom.score(theta, N, E)
negbinom.score(theta, N, E)
theta |
vector of parameter values. |
N |
vector of observed error profiles counts. |
E |
vector of expected error profiles counts. |
For further details see Myers et al. (2011).
negbinom.score
returns the vector of log-likelihood score values for the negative binomial distribution.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.
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.
dnegbinom
,
negbinom.em
,
negbinom.loglik
.
This function computes posterior predictive test statistics as described in Myers et al. (2011).
p.value(reps)
p.value(reps)
reps |
list of replications created with the |
For further details see Myers et al. (2011).
p-value
creates a list of p-values.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.constr.resamp
,
bhm.mcmc
,
bhm.resample
,
post.rep
.
## 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) reps <- post.rep(fit2, simdata) pvalues <- p.value(reps) ## End(Not run)
## 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) reps <- post.rep(fit2, simdata) pvalues <- p.value(reps) ## End(Not run)
Methods for function plot
in Package ‘graphics’ to be used with "mederrData", mederrFit and "mederrResample" objects.
signature(x = "mederrData", y = "missing")
Pictorial representation for a "mederrData" object.
signature(x = "mederrFit", y = "mederrFit")
Graphical representation of Markov Chain Monte Carlo simulations for a "mederrFit" object.
signature(x = "mederrResample", y = "missing")
Graphical representation of the resampling transformation for a "mederrResample" object.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
This function creates a list of replicated data for posterior predictive checking as described in Myers et al. (2011).
post.rep(model, dat)
post.rep(model, dat)
model |
|
dat |
an object of class "mederrData". |
For further details see Myers et al. (2011).
post.rep
returns a list of replicated data.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.constr.resamp
,
bhm.mcmc
,
bhm.resample
,
p.value
.
Random generation for a mixture of two negative binomial distributions.
rmixnegbinom(n, theta, E)
rmixnegbinom(n, theta, E)
n |
number of observations. |
theta |
vector of parameters for the negative binomial distribution mixture. |
E |
vector of (non-negative integer) expected counts. |
rmixnegbinom
generates random deviates corresponding to the E
and theta
values provided.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
Simulated data to use for illustrative purposes in the mederrRank
package.
data(simdata)
data(simdata)
An object of class mederrData
.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.
bhm.mcmc
,
mederrData
,
mederrFit
.
data("simdata", package = "mederrRank") summary(simdata) plot(simdata)
data("simdata", package = "mederrRank") summary(simdata) plot(simdata)
Methods for function summary
in Package ‘base’ to be used with "mederrData" and "mederrFit" objects.
signature(object = "mederrData")
Extracts summary information about the slots of a "mederrData" object.
signature(object = "mederrFit")
Extracts summary information about the slots of a "mederrFit" object.
Sergio Venturini [email protected],
Jessica A. Myers [email protected]
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.