R package prevalence | Center for Burden and Risk Assessment

Estimate true prevalence from individuals samples using multiple tests – conditional probability scheme

Description

Bayesian estimation of true prevalence from apparent prevalence obtained by applying multiple tests to individual samples. truePrevMulti implements the approach described by Berkvens et al. (2006), which uses a multinomial distribution to model observed test results, and in which conditional dependence between tests is modelled through conditional probabilities.

Usage

truePrevMulti(x, n, prior, nchains = 2, burnin = 10000, update = 10000,
              verbose = FALSE)

Arguments

`x`	Vector of apparent test results; see ‘Details’ below
`n`	The total sample size
`prior`	The prior distribution for `theta`; see ‘Details’ below
`nchains`	The number of chains used in the estimation process; 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

truePrevMulti calls on JAGS via the rjags package to estimate true prevalence from apparent prevalence in a Bayesian framework. truePrevMulti fits a multinomial model to the apparent test results obtained by testing individual samples with a given number of tests. To see the actual fitted model, see the model slot of the prev-class-object.

The vector of apparent tests results, x, must contain the number of samples corresponding to each combination of test results. To see how this vector is defined for the number of tests h at hand, use define_x.

The prior in the multinomial model consists of a vector theta, which holds values for the true prevalence (TP), the sensitivity and specificity of the first test (SE1, SP1), and the conditional dependencies between the results of the subsequent tests and the preceding one(s). To see how this vector is defined for the number of tests n at hand, use define_prior.

The values of prior can be specified in two ways, referred to as BUGS-style and list-style, respectively. See also below for some examples.

For BUGS-style specification, the values of theta should be given between curly brackets (i.e., {}), separated by line breaks. theta values can be specified to be deterministic (i.e., fixed), using the <- operator, or stochastic, using the ~ operator. In the latter case, the following distributions can be used:

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

Alternatively, theta values 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 parameterization.
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.

Value

A prev-class object.

Note

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

References

Berkvens D, Speybroeck N, Praet N, Adel A, Lesaffre E (2006) Estimating disease prevalence in a Bayesian framework using probabilistic constraints. Epidemiology 17:145-153. http://dx.doi.org/10.1097/01.ede.0000198422.64801.8d
Habib I, Sampers I, Uyttendaele M, De Zutter L, Berkvens D (2008) A Bayesian modelling framework to estimate Campylobacter prevalence and culture methods sensitivity: application to a chicken meat survey in Belgium. J Appl Microbiol 105:2002-2008. http://dx.doi.org/10.1111/j.1365-2672.2008.03902.x
Geurden T, Berkvens D, Casaert S, Vercruysse J, Claerebout E (2008) A Bayesian evaluation of three diagnostic assays for the detection of Giardia duodenalis in symptomatic and asymptomatic dogs. Vet Parasitol 157:14-20. http://dx.doi.org/10.1016/j.vetpar.2008.07.002

Examples

## ===================================================== ##
## 2-TEST EXAMPLE: Campylobacter                         ##
## ----------------------------------------------------- ##
## Two tests were performed on 656 chicken meat samples  ##
## -> T1 = enrichment culture                            ##
## -> T2 = direct plating                                ##
## The following assumption were made:                   ##
## -> TP is larger than 45% and smaller than 80%         ##
## -> SE1 must lie within 24% and 50%                    ##
## -> SP1 and SP2 both equal 100%                        ##
## -> beta(30, 12) describes P(T2+|D+,T1+)               ##
## The following results were obtained:                  ##
## -> 113 samples T1+,T2+                                ##
## ->  46 samples T1+,T2-                                ##
## -> 156 samples T1-,T2+                                ##
## -> 341 samples T1-,T2-                                ##
## ===================================================== ##

## how is the 2-test model defined?
define_x(2)
#> Definition of the apparent test results, 'x', for 2 tests:
#> x[1] : T1+,T2+ 
#> x[2] : T1+,T2- 
#> x[3] : T1-,T2+ 
#> x[4] : T1-,T2-

define_prior(2)
#> Conditional probability scheme
#> Definition of the prior, 'theta', for 2 tests: 
#> theta[1] : P(D+) = TP
#> theta[2] : P(T1+|D+) = SE1
#> theta[3] : P(T1-|D-) = SP1
#> theta[4] : P(T2+|D+,T1+)
#> theta[5] : P(T2+|D+,T1-)
#> theta[6] : P(T2-|D-,T1-)
#> theta[7] : P(T2-|D-,T1+)

## fit campylobacter 2-test model
campy <-
truePrevMulti(
  x = c(113, 46, 156, 341),
  n = 656,
  prior = {
    theta[1] ~ dunif(0.45, 0.80)
    theta[2] ~ dunif(0.24, 0.50)
    theta[3] <- 1
    theta[4] ~ dbeta(30, 12)
    theta[5] ~ dbeta(1, 1)
    theta[6] <- 1
    theta[7] <- 1
  }
)

## fit same model using 'list-style'
campy <-
truePrevMulti(
  x = c(113, 46, 156, 341),
  n = 656,
  prior =
    list(
      theta1 = list(dist = "uniform", min = 0.45, max = 0.80),
      theta2 = list(dist = "uniform", min = 0.24, max = 0.50),
      theta3 = 1,
      theta4 = list(dist = "beta", alpha = 30, beta = 12),
      theta5 = list(dist = "beta", alpha = 1, beta = 1),
      theta6 = 1,
      theta7 = 1
    )
)

## show model results
campy
#>      mean median  mode    sd  2.5% 97.5%
#>  TP 0.618  0.607 0.538 0.090 0.480 0.787
#> SE1 0.398  0.398 0.362 0.058 0.295 0.493
#> SP1 1.000  1.000 1.000 0.000 1.000 1.000
#> SE2 0.677  0.675 0.609 0.099 0.509 0.847
#> SP2 1.000  1.000 1.000 0.000 1.000 1.000
#> 
#> Multivariate BGR statistic = 1.001
#> BGR values substantially above 1 indicate lack of convergence
#> Bayes-P statistic = 0.47 
#> Bayes-P values substantially different from 0.5 indicate lack of convergence

## explore model structure
str(campy)         # overall structure
#> Formal class 'prev' [package "prevalence"] with 4 slots
#>   ..@ par        :List of 7
#>   .. ..$ x      : num [1:4] 113 46 156 341
#>   .. ..$ n      : num 656
#>   .. ..$ prior  :List of 7
#>   .. .. ..$ :List of 2
#>   .. .. .. ..$ d: chr "uniform"
#>   .. .. .. ..$ p: num [1:2] 0.45 0.8
#>   .. .. ..$ :List of 2
#>   .. .. .. ..$ d: chr "uniform"
#>   .. .. .. ..$ p: num [1:2] 0.24 0.5
#>   .. .. ..$ :List of 2
#>   .. .. .. ..$ d: chr "fixed"
#>   .. .. .. ..$ p: num 1
#>   .. .. ..$ :List of 2
#>   .. .. .. ..$ d: chr "beta"
#>   .. .. .. ..$ p: num [1:2] 30 12
#>   .. .. ..$ :List of 2
#>   .. .. .. ..$ d: chr "beta"
#>   .. .. .. ..$ p: num [1:2] 1 1
#>   .. .. ..$ :List of 2
#>   .. .. .. ..$ d: chr "fixed"
#>   .. .. .. ..$ p: num 1
#>   .. .. ..$ :List of 2
#>   .. .. .. ..$ d: chr "fixed"
#>   .. .. .. ..$ p: num 1
#>   .. ..$ nchains: num 2
#>   .. ..$ burnin : num 10000
#>   .. ..$ update : num 10000
#>   .. ..$ inits  : NULL
#>   ..@ model      :Class 'prevModel'  chr [1:31] "model {" "x[1:4] ~ dmulti(AP[1:4], n)" "" "AP[1] <- theta[1]*theta[2]*theta[4] + (1-theta[1])*(1-theta[3])*(1-theta[7])" ...
#>   ..@ mcmc       :List of 6
#>   .. ..$ TP    :List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.632 0.602 0.605 0.618 0.549 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.637 0.643 0.696 0.737 0.735 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   .. ..$ SE1   :List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.399 0.406 0.405 0.362 0.409 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.377 0.399 0.388 0.306 0.301 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   .. ..$ SP1   :List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   .. ..$ SE2   :List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.66 0.649 0.729 0.646 0.674 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.659 0.649 0.626 0.595 0.5 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   .. ..$ SP2   :List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   .. ..$ bayesP:List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 1 0 1 0 1 0 0 0 0 1 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 1 0 1 1 1 1 1 1 0 0 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   ..@ diagnostics:List of 3
#>   .. ..$ DIC   :List of 3
#>   .. .. ..$ deviance:Class 'mcarray'  Named num 21.2
#>   .. .. .. .. ..- attr(*, "names")= chr "x"
#>   .. .. ..$ penalty :Class 'mcarray'  Named num 2.67
#>   .. .. .. .. ..- attr(*, "names")= chr "x"
#>   .. .. ..$ type    : chr "pD"
#>   .. .. ..- attr(*, "class")= chr "dic"
#>   .. ..$ BGR   :List of 2
#>   .. .. ..$ psrf : num [1:3, 1:2] 1 1 1 1 1 ...
#>   .. .. .. ..- attr(*, "dimnames")=List of 2
#>   .. .. .. .. ..$ : chr [1:3] "SE1" "SE2" "TP"
#>   .. .. .. .. ..$ : chr [1:2] "Point est." "Upper C.I."
#>   .. .. ..$ mpsrf: num 1
#>   .. .. ..- attr(*, "class")= chr "gelman.diag"
#>   .. ..$ bayesP: num 0.466
str(campy@par)     # structure of slot 'par'
#> List of 7
#>  $ x      : num [1:4] 113 46 156 341
#>  $ n      : num 656
#>  $ prior  :List of 7
#>   ..$ :List of 2
#>   .. ..$ d: chr "uniform"
#>   .. ..$ p: num [1:2] 0.45 0.8
#>   ..$ :List of 2
#>   .. ..$ d: chr "uniform"
#>   .. ..$ p: num [1:2] 0.24 0.5
#>   ..$ :List of 2
#>   .. ..$ d: chr "fixed"
#>   .. ..$ p: num 1
#>   ..$ :List of 2
#>   .. ..$ d: chr "beta"
#>   .. ..$ p: num [1:2] 30 12
#>   ..$ :List of 2
#>   .. ..$ d: chr "beta"
#>   .. ..$ p: num [1:2] 1 1
#>   ..$ :List of 2
#>   .. ..$ d: chr "fixed"
#>   .. ..$ p: num 1
#>   ..$ :List of 2
#>   .. ..$ d: chr "fixed"
#>   .. ..$ p: num 1
#>  $ nchains: num 2
#>  $ burnin : num 10000
#>  $ update : num 10000
#>  $ inits  : NULL
str(campy@mcmc)    # structure of slot 'mcmc'
#> List of 6
#>  $ TP    :List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.632 0.602 0.605 0.618 0.549 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.637 0.643 0.696 0.737 0.735 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
#>  $ SE1   :List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.399 0.406 0.405 0.362 0.409 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.377 0.399 0.388 0.306 0.301 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
#>  $ SP1   :List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
#>  $ SE2   :List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.66 0.649 0.729 0.646 0.674 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.659 0.649 0.626 0.595 0.5 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
#>  $ SP2   :List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 1 1 1 1 1 1 1 1 1 1 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
#>  $ bayesP:List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 1 0 1 0 1 0 0 0 0 1 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 1 0 1 1 1 1 1 1 0 0 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
campy@model        # fitted model
#> model {
#>   x[1:4] ~ dmulti(AP[1:4], n)
#>   
#>   AP[1] <- theta[1]*theta[2]*theta[4] + (1-theta[1])*(1-theta[3])*(1-theta[7])
#>   AP[2] <- theta[1]*theta[2]*(1-theta[4]) + (1-theta[1])*(1-theta[3])*theta[7]
#>   AP[3] <- theta[1]*(1-theta[2])*theta[5] + (1-theta[1])*theta[3]*(1-theta[6])
#>   AP[4] <- theta[1]*(1-theta[2])*(1-theta[5]) + (1-theta[1])*theta[3]*theta[6]
#>   
#>   theta[1] ~ dunif(0.45, 0.8)
#>   theta[2] ~ dunif(0.24, 0.5)
#>   theta[3] <- 1
#>   theta[4] ~ dbeta(30, 12)
#>   theta[5] ~ dbeta(1, 1)
#>   theta[6] <- 1
#>   theta[7] <- 1
#>   
#>   x2[1:4] ~ dmulti(AP[1:4], n)
#>   for (i in 1:4) {
#>     d1[i] <- x[i] * log(max(x[i],1) / (AP[i]*n))
#>     d2[i] <- x2[i] * log(max(x2[i],1) / (AP[i]*n))
#>   }
#>   G0 <- 2 * sum(d1[])
#>   Gt <- 2 * sum(d2[])
#>   bayesP <- step(G0 - Gt)
#>   
#>   TP <- theta[1]
#>   SE1 <- theta[2]
#>   SP1 <- theta[3]
#>   SE2 <- theta[2] * theta[4] + (1-theta[2]) * theta[5]
#>   SP2 <- theta[3] * theta[6] + (1-theta[3]) * theta[7]
#> }
campy@diagnostics  # DIC, BGR and Bayes-P values
#> $DIC
#> Mean deviance:  21.17 
#> penalty 2.672 
#> Penalized deviance: 23.84 
#> 
#> $BGR
#> Potential scale reduction factors:
#> 
#>     Point est. Upper C.I.
#> SE1          1          1
#> SE2          1          1
#> TP           1          1
#> 
#> Multivariate psrf
#> 
#> 1
#> 
#> $bayesP
#> [1] 0.46595

## standard methods
print(campy)
#>      mean median  mode    sd  2.5% 97.5%
#>  TP 0.618  0.607 0.538 0.090 0.480 0.787
#> SE1 0.398  0.398 0.362 0.058 0.295 0.493
#> SP1 1.000  1.000 1.000 0.000 1.000 1.000
#> SE2 0.677  0.675 0.609 0.099 0.509 0.847
#> SP2 1.000  1.000 1.000 0.000 1.000 1.000
#> 
#> Multivariate BGR statistic = 1.001
#> BGR values substantially above 1 indicate lack of convergence
#> Bayes-P statistic = 0.47 
#> Bayes-P values substantially different from 0.5 indicate lack of convergence
summary(campy)
#> $TP
#>                 mean    median      mode         sd         var      2.5%
#> chain 1    0.6199229 0.6078247 0.5368791 0.09150965 0.008374015 0.4812390
#> chain 2    0.6158595 0.6061822 0.5360000 0.08914862 0.007947477 0.4787141
#> all chains 0.6178912 0.6070227 0.5375477 0.09035743 0.008164466 0.4802267
#>                97.5% samples
#> chain 1    0.7873944   10000
#> chain 2    0.7867800   10000
#> all chains 0.7871350   20000
#> 
#> $SE1
#>                 mean    median      mode         sd         var      2.5%
#> chain 1    0.3969175 0.3974266 0.4507818 0.05872172 0.003448240 0.2934889
#> chain 2    0.3986561 0.3993997 0.3954451 0.05691255 0.003239038 0.2954941
#> all chains 0.3977868 0.3984066 0.3617115 0.05782930 0.003344227 0.2945072
#>                97.5% samples
#> chain 1    0.4932038   10000
#> chain 2    0.4931574   10000
#> all chains 0.4931970   20000
#> 
#> $SP1
#>            mean median mode sd var 2.5% 97.5% samples
#> chain 1       1      1    1  0   0    1     1   10000
#> chain 2       1      1    1  0   0    1     1   10000
#> all chains    1      1    1  0   0    1     1   20000
#> 
#> $SE2
#>                 mean    median      mode         sd         var      2.5%
#> chain 1    0.6749214 0.6738424 0.6107965 0.09985975 0.009971970 0.5092623
#> chain 2    0.6781265 0.6771043 0.6639961 0.09719916 0.009447677 0.5086821
#> all chains 0.6765239 0.6754320 0.6089149 0.09854900 0.009711906 0.5090104
#>                97.5% samples
#> chain 1    0.8473936   10000
#> chain 2    0.8461685   10000
#> all chains 0.8468951   20000
#> 
#> $SP2
#>            mean median mode sd var 2.5% 97.5% samples
#> chain 1       1      1    1  0   0    1     1   10000
#> chain 2       1      1    1  0   0    1     1   10000
#> all chains    1      1    1  0   0    1     1   20000
par(mfrow = c(2, 2))
plot(campy)         # shows plots of TP by default

plot(campy, "SE1")  # same plots for SE1

plot(campy, "SE2")  # same plots for SE2

## coda plots of TP, SE1, SE2
par(mfrow = c(1, 3))
densplot(campy, col = "red")

traceplot(campy)

gelman.plot(campy)

autocorr.plot(campy)

## ===================================================== ##
## 3-TEST EXAMPLE: Giardia                               ##
## ----------------------------------------------------- ##
## Three tests were performed on stools from 272 dogs    ##
## -> T1 = immunofluorescence assay                      ##
## -> T2 = direct microscopy                             ##
## -> T3 = SNAP immunochromatography                     ##
## The following assumption were made:                   ##
## -> TP is smaller than 20%                             ##
## -> SE1 must be higher than 80%                        ##
## -> SP1 must be higher than 90%                        ##
## The following results were obtained:                  ##
## ->   6 samples T1+,T2+,T3+                            ##
## ->   4 samples T1+,T2+,T3-                            ##
## ->  12 samples T1+,T2-,T3+                            ##
## ->  12 samples T1+,T2-,T3-                            ##
## ->   1 sample  T1-,T2+,T3+                            ##
## ->  14 samples T1-,T2+,T3-                            ##
## ->   3 samples T1-,T2-,T3+                            ##
## -> 220 samples T1-,T2-,T3-                            ##
## ===================================================== ##

## how is the 3-test model defined?
define_x(3)
#> Definition of the apparent test results, 'x', for 3 tests:
#> x[1] : T1+,T2+,T3+ 
#> x[2] : T1+,T2+,T3- 
#> x[3] : T1+,T2-,T3+ 
#> x[4] : T1+,T2-,T3- 
#> x[5] : T1-,T2+,T3+ 
#> x[6] : T1-,T2+,T3- 
#> x[7] : T1-,T2-,T3+ 
#> x[8] : T1-,T2-,T3-

define_prior(3)
#> Conditional probability scheme
#> Definition of the prior, 'theta', for 3 tests: 
#> theta[1] : P(D+) = TP
#> theta[2] : P(T1+|D+) = SE1
#> theta[3] : P(T1-|D-) = SP1
#> theta[4] : P(T2+|D+,T1+)
#> theta[5] : P(T2+|D+,T1-)
#> theta[6] : P(T2-|D-,T1-)
#> theta[7] : P(T2-|D-,T1+)
#> theta[8] : P(T3+|D+,T1+,T2+)
#> theta[9] : P(T3+|D+,T1+,T2-)
#> theta[10] : P(T3+|D+,T1-,T2+)
#> theta[11] : P(T3+|D+,T1-,T2-)
#> theta[12] : P(T3-|D-,T1-,T2-)
#> theta[13] : P(T3-|D-,T1-,T2+)
#> theta[14] : P(T3-|D-,T1+,T2-)
#> theta[15] : P(T3-|D-,T1+,T2+)

## fit giardia 3-test model
giardia <-
truePrevMulti(
  x = c(6, 4, 12, 12, 1, 14, 3, 220),
  n = 272,
  prior = {
    theta[1] ~ dunif(0.00, 0.20)
    theta[2] ~ dunif(0.90, 1.00)
    theta[3] ~ dunif(0.80, 1.00)
    theta[4] ~ dbeta(1, 1)
    theta[5] ~ dbeta(1, 1)
    theta[6] ~ dbeta(1, 1)
    theta[7] ~ dbeta(1, 1)
    theta[8] ~ dbeta(1, 1)
    theta[9] ~ dbeta(1, 1)
    theta[10] ~ dbeta(1, 1)
    theta[11] ~ dbeta(1, 1)
    theta[12] ~ dbeta(1, 1)
    theta[13] ~ dbeta(1, 1)
    theta[14] ~ dbeta(1, 1)
    theta[15] ~ dbeta(1, 1)
  }
)

## show model results
giardia
#>      mean median  mode    sd  2.5% 97.5%
#>  TP 0.069  0.068 0.055 0.037 0.006 0.144
#> SE1 0.950  0.950 0.942 0.029 0.903 0.997
#> SP1 0.931  0.932 0.933 0.036 0.862 0.994
#> SE2 0.361  0.321 0.264 0.212 0.051 0.886
#> SP2 0.917  0.918 0.920 0.020 0.875 0.953
#> SE3 0.521  0.524 0.527 0.173 0.160 0.860
#> SP3 0.943  0.944 0.948 0.024 0.894 0.983
#> 
#> Multivariate BGR statistic = 1.0014
#> BGR values substantially above 1 indicate lack of convergence
#> Bayes-P statistic = 0.46 
#> Bayes-P values substantially different from 0.5 indicate lack of convergence

## coda densplots
par(mfcol = c(2, 4))
densplot(giardia, col = "red")