rrisk.BayesPEM: Bayesian Prevalence Estimation under Misclassification (PEM)

In BfRstats/rriskBayes2: Predefined Bayes models fitted with Markov chain Monte Carlo (MCMC) (related to the 'rrisk' project)

Description

Bayesian PEM models provide the posterior distribution for the true prevalence (pi), diagnostic sensitivity (se) and specificity (sp) for a given empirical prevalence estimate using physically pooled samples (if k>1) and priors for the model parameters. The misclassification parameters (se and sp) can be specified at the level of the pool or individual level of testing. On the other side, the function estimates the true prevalence based on the results (x/n) of an application study with individual samples (if k=1) using a diagnostic test, for which some prior information on sensitivity and/or specificity is available.

Usage

rrisk.BayesPEM(x, n, k, simulation = FALSE, 
 prior.pi = c(1, 1), prior.se, prior.sp,
 misclass = "pool", chains = 3, burn = 4000, 
 thin = 1, update = 5000, 
 max.time = "15minutes", plots = FALSE)

Arguments

x	Scalar value for number of pools (k>1) or individual outcomes (k=1) with positive test result.
n	Scalar value for number of pools tested (k>1) or the sample size in application study (k=1).
k	Scalar value for number of individual samples physically combined into one pool; set k>1 for pooled sampling and k=1 for individual sampling (default 1).
prior.pi	Numeric vector containing parameters of a beta distribution as prior for prevalence pi, e.g., pi ~ prior.pi(*,*) = beta(*,*).
prior.se	Numeric vector containing parameters of a beta distribution as prior for sensitivity se, e.g., se ~ prior.se(*,*) = beta(*,*). For fixed sensitivity scalar value.
prior.sp	Numeric vector containing parameters of a beta distribution as prior for specificity sp, e.g., sp ~ prior.sp(*,*) = beta(*,*). For fixed specifity scalar value.
simulation	Not used any longer.
misclass	Character with legal character entries: individual, individual-fix-se, individual-fix-sp, individual-fix-se-sp, pool, pool-fix-se, pool-fix-sp or pool-fix-se-sp. fix-se: fixed sensitivity fix-sp: fixed specifity fix-se-sp: fixed sensitivity AND fixed specifity.
chains	Positive single numeric value, number of independent MCMC chains (default 3).
burn	Positive single numeric value, length of the burn-in period (default 4000).
thin	Positive single numeric value (default 1). The samples from every k-th iteration will be used for inference, where k is the value of thin. Setting thin > 1 can help to reduce the autocorrelation in the sample.
update	Positive single numeric value, length of update iterations for estimation (default 5000).
max.time	Maximum time for which the function is allowed to extend the chains. Acceptable units include 'seconds', 'minutes', 'hours', 'days', 'weeks' (default "15minutes") (see autorun.jags).
plots	Logical, if TRUE the diagnostic plots will be displayed.

Details

The application data (k=1) has one degree of freedom while the underlying model has three unknown parameters. Thus, the model is not identifiable and informative priors on at least two model parameters are required. The Bayesian models for estimation prevalence, sensitivity and specificity take the following forms in rjags/JAGS (originally BRugs/Winbugs) syntax:

Misclassifications at the individual level (k=1 and misclass="individual")

model{

      pi  ~ dbeta(prior.pi[1], prior.pi[2])

      se ~ dbeta(prior.se[1], prior.se[2])

      sp ~ dbeta(prior.sp[1], prior.sp[2])

      ap <- pi * se + (1-pi) * (1-sp)
      
      x ~ dbin(ap, n)
   }

Misclassifications at the individual level with fixed sensitivity resp. specifity (k=1 and misclass="individual-fix-sp") Both se and sp could be set to fixed values.

model{
      pi ~ dbeta(prior.pi[1], prior.pi[2])

      se resp. sp ~ dbeta(prior.sp[1], prior.sp[2])

      sp resp. se <-  fix

      ap <- pi*se + (1-pi)*(1-sp)

      x ~ dbin(p,n)
   }

For misclassification at the pool-level (k>1 and misclass="pool")

model{

      pi ~ dbeta(prior.pi[1], prior.pi[2])

      se ~ dbeta(prior.se[1], prior.se[2])

      sp ~ dbeta(prior.sp[1], prior.sp[2])

      p.neg <- pow(1-pi, k)

      p <- (1-p.neg)*se + p.neg*(1-sp)

      x ~ dbin(p, n)
     }

For misclassification at the pool-level (k>1) and fixed sensitivity and specifity (misclass="pool-fix-se-sp")

model{

      pi ~ dbeta(prior.pi[1], prior.pi[2])

      se <- fix

      sp <- fix

      p.neg <- pow(1-pi,k)

      p <- (1-p.neg)*se + p.neg*(1-sp)

      x ~ dbin(p,n)
     }

and for comparison (k>1)

model{

      pi1 ~ dbeta(prior.pi[1], prior.pi[2])

      pi2 ~ dbeta(prior.pi[1], prior.pi[2])

      se ~ dbeta(prior.se[1], prior.se[2])

      sp ~ dbeta(prior.sp[1], prior.sp[2])

      x1 <- x

      x2 <- x

      p.neg <- pow(1-pi1, k)

      p.pos <- (1-p.neg)*se + p.neg*(1-sp)

      x1 ~ dbin(p.pos, n)

      ap <- pi2*se + (1-pi2)*(1-sp)

      p <- 1- pow(1-ap,k)

      x2 ~ dbin(p,n)
     }

Value

The function rrisk.BayesZIP returns an instance of the bayesmodelClass class containing following information:

convergence	Logical, whether the model has converged (assessed by the user).
results	Data frame containing statitsics of the posterior distribution.
jointpost	Data frame giving the joint posterior probability distribution.
nodes	Names of the parameters jointly estimated by the Bayes model.
model	Model in rjags/JAGS syntax as a character string.
chains	Number of independent MCMC chains.
burn	Length of burn-in period.
update	Length of update iterations for estimation.

References

Greiner, M., Belgoroski, N., NN (2011). Estimating prevalence using diagnostic data with pooled samples: R function rrisk.BayesPEM. J.Stat.Software (in preparation).

Cowling, D.W., Gardner, I.A. and Johnson W.O. (1999). Comparison of methods for estimation of individual-level prevalence based on pooled samples, Prev.Vet.Med. 39: 211-225.

Rogan, W.J. and B. Gladen (1978). Estimating prevalence from the results of a screening test. Am. J. Epidemiol. 107: 71-76.

Examples

------------------------------------------
Example of PEM model
------------------------------------------
# generate PEM data at individual level

n <- 100
x <- 14
k <- 1
pi_prior <- c(1, 1)
se_prior <- c(64, 4)
sp_prior <- c(94, 16)
misclass <- "individual"

# run model

resPEM <- rrisk.BayesPEM(x = x, n = n, k = k, 
prior.pi = pi_prior, prior.se = se_prior, prior.sp = sp_prior,  
misclass = misclass, plots=TRUE)
  
# generate data for pooled sampling and fixed sensitivity and specifity

n <- 100
x <- 14
k <- 4
pi_prior <- c(1, 1)
sp_prior <- 1
se_prior <- 0.7
misclass <- "pool-fix-se-sp"

# run model

resPEM <- rrisk.BayesPEM(x = x, n = n, k = k, 
prior.pi = pi_prior, prior.se = se_prior, prior.sp = sp_prior, 
misclass = misclass, plots=TRUE)

BfRstats/rriskBayes2 documentation built on May 5, 2019, 2:42 p.m.