R package prevalence | Center for Burden and Risk Assessment

Estimate true prevalence from pooled samples

Description

Bayesian estimation of true prevalence from apparent prevalence obtained by testing pooled samples.

Usage

truePrevPools(x, n, SE = 1, SP = 1, prior = c(1, 1),
              nchains = 2, burnin = 10000, update = 10000,
              verbose = FALSE)

Arguments

`x`	The vector of indicator variables, indicating whether a pool was positive (`"1"`) or negative (`"0"`)
`n`	The vector of pool sizes
`SE, SP`	The prior distribution for sensitivity (SE) and specificity (SP); see ‘Details’ below for specification of these distributions
`prior`	The parameters of the prior Beta distribution for true prevalence; defaults to `c(1, 1)`
`nchains`	The number of chains used in the estimation process; `nchains` must be 2
`burnin`	The number of discarded model iterations; defaults to 10,000
`update`	The number of withheld model iterations; defaults to 10,000
`verbose`	Logical flag, indicating if JAGS process output should be printed to the R console; defaults to `FALSE`

Details

truePrevPools calls on JAGS via rjags to estimate the true prevalence from the apparent prevalence in a Bayesian framework. The default model, in BUGS language, is given below. To see the actual fitted model, see the model slot of the prev-class object.

model {
    for (i in 1:N) {
      x[i] ~ dbern(AP[i])
      AP[i] <- SEpool[i] * (1 - pow(1 - TP, n[i])) + (1 - SPpool[i]) * pow(1 - TP, n[i])
      SEpool[i] <- 1 - (pow(1 - SE, n[i] * TP) * pow(SP, n[i] * (1 - TP)))
      SPpool[i] <- pow(SP, n[i])
    }
  # SE ~ user-defined
  # SP ~ user-defined
  TP ~ dbeta(prior[1], prior[2])
  }

The test sensitivity (SE) and specificity (SP) can be specified by the user, independently, as one of "fixed" (default), "uniform", "beta", "pert", or "beta-expert". Note that SE and SP must correspond to the test characteristics for testing individual samples; truePrevPools will calculate SEpool and SPpool, the sensitivity and specificitiy for testing pooled samples, based on Boelaert et al. (2000).

Distribution parameters can be specified in a named list() as follows:

Fixed: list(dist = "fixed", par)
Uniform: list(dist = "uniform", min, max)
Beta: list(dist = "beta", alpha, beta)
PERT: list(dist = "pert", method, a, m, b, k)
'method' must be "classic" or "vose";
'a' denotes the pessimistic (minimum) estimate, 'm' the most likely estimate, and 'b' the optimistic (maximum) estimate;
'k' denotes the scale parameter.
See betaPERT for more information on Beta-PERT parametrization.
Beta-Expert: list(dist = "beta-expert", mode, mean, lower, upper, p)
'mode' denotes the most likely estimate, 'mean' the mean estimate;
'lower' denotes the lower bound, 'upper' the upper bound;
'p' denotes the confidence level of the expert.
Only mode or mean should be specified; lower and upper can be specified together or alone.
See betaExpert for more information on Beta-Expert parameterization.

For Uniform, Beta and Beta-PERT distributions, BUGS-style short-hand notation is also allowed:

Uniform: ~dunif(min, max)
Beta: ~dbeta(alpha, beta)
Beta-PERT: ~dpert(min, mode, max)

Value

An object of class prev-class.

Note

Markov chain Monte Carlo sampling in truePrev is performed by JAGS (Just Another Gibbs Sampler) through the rjags package. JAGS can be downloaded from http://sourceforge.net/projects/mcmc-jags/.

References

Speybroeck N, Williams CJ, Lafia KB, Devleesschauwer B, Berkvens D (2012) Estimating the prevalence of infections in vector populations using pools of samples. Med Vet Entomol 26:361-371. http://dx.doi.org/10.1111/j.1365-2915.2012.01015.x
Boelaert F, Walravens K, Biront P, Vermeersch JP, Berkvens D, Godfroid J (2000) Prevalence of paratuberculosis (Johne’s disease) in the Belgian cattle population. Vet Microbiol 77:269-281. http://dx.doi.org/10.1016/S0378-1135(00)00312-6

Examples

## Sandflies in Aurabani, Nepal, 2007
pool_results <- c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0)
pool_sizes <- c(2, 1, 6, 10, 1, 7, 1, 4, 1, 3)

## Sensitivity ranges uniformly between 60% and 95%
## Specificity is considered to be 100%

## .. BUGS-style
truePrevPools(x = pool_results, n = pool_sizes,
              SE = ~dunif(0.60, 0.95), SP = 1)
#>     mean median  mode    sd  2.5% 97.5%
#> TP 0.100  0.093 0.083 0.050 0.025 0.218
#> SE 0.764  0.759 0.630 0.101 0.607 0.941
#> SP 1.000  1.000 1.000 0.000 1.000 1.000
#> 
#> Multivariate BGR statistic = 1
#> BGR values substantially above 1 indicate lack of convergence

## .. list-style
SE <- list(dist = "uniform", min = 0.60, max = 0.95)
truePrevPools(x = pool_results, n = pool_sizes,
              SE = SE, SP = 1)
#>     mean median  mode    sd  2.5% 97.5%
#> TP 0.099  0.092 0.079 0.050 0.025 0.217
#> SE 0.765  0.760 0.633 0.101 0.607 0.940
#> SP 1.000  1.000 1.000 0.000 1.000 1.000
#> 
#> Multivariate BGR statistic = 1.0001
#> BGR values substantially above 1 indicate lack of convergence