Estimate true prevalence from individuals samples

Share:

Description

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

Usage

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

Arguments

x

The apparent number of positive samples

n

The sample size

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; 'n' 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

truePrev calls on JAGS/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-object.

1
2
3
4
5
6
7
8
  model {
    x ~ dbin(AP, n)
    AP <- SE * TP + (1 - SP) * (1 - TP)
    # SE ~ user-defined (see below)
    # SP ~ user-defined (see below)
    TP ~ dbeta(prior[1], prior[2])
  }
  

The test sensitivity (SE) and specificity (SP) can be specified, independently, as one of "fixed", "uniform", "beta", "pert", or "beta-expert", with "fixed" as the default.

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)

  • 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 parametrization.

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.

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/.

Author(s)

Brecht Devleesschauwer <brechtdv@gmail.com>

References

See Also

coda for various functions that can be applied to the prev@mcmc object
truePrevMulti: estimate true prevalence from apparent prevalence obtained by testing individual samples with multiple tests, using a conditional probability scheme
truePrevMulti2: estimate true prevalence from apparent prevalence obtained by testing individual samples with multiple tests, using a covariance scheme
truePrevPools: estimate true prevalence from apparent prevalence obtained by testing pooled samples
betaPERT: calculate the parameters of a Beta-PERT distribution
betaExpert: calculate the parameters of a Beta distribution based on expert opinion

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
## Taenia solium cysticercosis in Nepal
## 142 positives out of 742 pigs sampled

## Model SE and SP based on literature data
## Sensitivity ranges uniformly between 60% and 100%
## Specificity ranges uniformly between 75% and 100%
#> BUGS-style:
truePrev(x = 142, n = 742,
         SE = ~dunif(0.60, 1.00), SP = ~dunif(0.75, 1.00))

#> list-style:
SE <- list(dist = "uniform", min = 0.60, max = 1.00)
SP <- list(dist = "uniform", min = 0.75, max = 1.00)
truePrev(x = 142, n = 742, SE = SE, SP = SP)

## Model SE and SP based on expert opinions
## Sensitivity lies in between 60% and 100%; most likely value is 90%
## Specificity is with 95% confidence larger than 75%; most likely value is 90%
SE <- list(dist = "pert", a = 0.60, m = 0.90, b = 1.00)
SP <- list(dist = "beta-expert", mode = 0.90, lower = 0.75, p = 0.95)
truePrev(x = 142, n = 742, SE = SE, SP = SP)

## Model SE and SP as fixed values (each 90%)
truePrev(x = 142, n = 742, SE = 0.90, SP = 0.90)