the prevalence package
tools for prevalence assessment studies.

‘prev’ class structure and methods


Functions truePrev, truePrevPools, truePrevMulti and truePrevMulti2 create objects of S4 class prev. Different methods are available for manipulating these objects.

Slots

Objects of class prev are so-called S4 objects. You can think of these objects as being formally defined lists. Instead of using the $-operator for accessing the different elements, which are called slots, the @-operator needs to be used. Objects of class prev contain the following four slots:

@par A list of input parameters.
@model The fitted Bayesian model, in BUGS language (and given S3 class prevModel).
@mcmc A list, with one element per chain, of the simulated true prevalences (and sensitivities and specificities in the case of truePrevMulti and truePrevMulti2).
@diagnostics A list with elements for the Deviance Information Criterion ($DIC), the Brooks-Gelman-Rubin statistic ($BGR), and in the case of truePrevMulti and truePrevMulti2, the Bayes-P statistic ($bayesP)

General methods

print(x, conf.level = 0.95, dig = 3, ...)
x An object of class prev
conf.level Confidence level to be used in credibility interval
dig Number of decimal digits to print
... Other arguments to pass to the print function

The print method prints the mean, median, mode, standard deviation and credibility interval of the estimated true prevalence (and sensitivities and specificities, for prev objects created by truePrevMulti and truePrevMulti2). In addition, the Brooks-Gelman-Rubin statistic and corresponding upper confidence limit are printed, or the multivariate BGR statistic for prev objects created by truePrevMulti and truePrevMulti2. BGR values substantially above 1 indicate lack of convergence. For prev objects created by truePrevMulti, the Bayes-P statistic is also printed. Bayes-P should be as close to 0.5 as possible.

summary method

summary(object, conf.level)
object An object of class prev
conf.level Confidence level to be used in credibility interval

The summary method returns the mean, median, mode, standard deviation, variance, credibility interval and number of samples for each chain separately and for all chains combined.

as.matrix method

as.matrix(x, iters = FALSE, chains = FALSE)
x An object of class prev
iters Logical flag, indicating whether a column should be added for iteration number; defaults to FALSE
chains Logical flag, indicating whether a column should be added for chain number; defaults to FALSE

The as.matrix method converts objects of class prev to matrix.

plot method

plot(x, y = NULL, ...)
x An object of class prev
y Which parameter to plot? Defaults to NULL
... Other arguments to pass to the plot function

The plot method generates density, trace, Brooks-Gelman-Rubin and autocorrelation plots.

Methods inherited from the coda package

densplot method

densplot(x, exclude_fixed = TRUE, ...)
x An object of class prev
exclude_fixed Should fixed parameters be excluded from plotting? Defaults to TRUE
... Other arguments to pass to the densplot function

Displays a plot of the density estimate for each variable in x, calculated by the density function.

traceplot method

traceplot(x, exclude_fixed = TRUE, ...)
x An object of class prev
exclude_fixed Should fixed parameters be excluded from plotting? Defaults to TRUE
... Other arguments to pass to the traceplot function

Displays a plot of iterations versus sampled values for each variable in the chain, with a separate plot per variable.

autocorr.plot method

autocorr.plot(x, exclude_fixed = TRUE, chain = 1, ...)
x An object of class prev
exclude_fixed Should fixed parameters be excluded from plotting? Defaults to TRUE
chain Which chain to plot in autocorr.plot? Defaults to 1
... Other arguments to pass to the autocorr.plot function

Plots the autocorrelation function for each variable in each chain in x.

Examples

## Taenia solium cysticercosis in Nepal
SE <- list(dist = "uniform", min = 0.60, max = 1.00)
SP <- list(dist = "uniform", min = 0.75, max = 1.00)
TP <- truePrev(x = 142, n = 742, SE = SE, SP = SP)
## Summarize estimates per chain
summary(TP)
#> $TP
#>                 mean    median      mode         sd         var
#> chain 1    0.1371857 0.1371797 0.1899004 0.07788892 0.006066685
#> chain 2    0.1336157 0.1321198 0.1831823 0.07938643 0.006302206
#> all chains 0.1354007 0.1344100 0.1881068 0.07865953 0.006187322
#>                   2.5%     97.5% samples
#> chain 1    0.007131391 0.2847301   10000
#> chain 2    0.006826268 0.2867341   10000
#> all chains 0.006953099 0.2858581   20000
#> 
#> $SE
#>                 mean    median      mode        sd        var      2.5%
#> chain 1    0.7781083 0.7664154 0.6305376 0.1160362 0.01346441 0.6066010
#> chain 2    0.7806555 0.7695690 0.6425182 0.1150805 0.01324352 0.6082347
#> all chains 0.7793819 0.7678970 0.6302666 0.1155635 0.01335492 0.6071912
#>                97.5% samples
#> chain 1    0.9867098   10000
#> chain 2    0.9874131   10000
#> all chains 0.9870292   20000
#> 
#> $SP
#>                 mean    median      mode         sd         var      2.5%
#> chain 1    0.9009873 0.9005704 0.8490157 0.05666717 0.003211168 0.8045467
#> chain 2    0.8990168 0.8968261 0.8380739 0.05781682 0.003342784 0.8032723
#> all chains 0.9000021 0.8986753 0.8478642 0.05725192 0.003277783 0.8039108
#>                97.5% samples
#> chain 1    0.9938295   10000
#> chain 2    0.9939376   10000
#> all chains 0.9938842   20000
## Diagnostic plots
par(mfrow = c(2, 2))
plot(TP)

## Generic plots from package coda
par(mfrow = c(1, 3))
densplot(TP)


par(mfrow = c(1, 3))
traceplot(TP)


par(mfrow = c(1, 3))
gelman.plot(TP)


par(mfrow = c(1, 3))
autocorr.plot(TP)

## Use 'str()' to see the structure of object TP
str(TP)
#> Formal class 'prev' [package "prevalence"] with 4 slots
#>   ..@ par        :List of 9
#>   .. ..$ x      : num 142
#>   .. ..$ n      : num 742
#>   .. ..$ SE     :List of 2
#>   .. .. ..$ d: chr "uniform"
#>   .. .. ..$ p: num [1:2] 0.6 1
#>   .. ..$ SP     :List of 2
#>   .. .. ..$ d: chr "uniform"
#>   .. .. ..$ p: num [1:2] 0.75 1
#>   .. ..$ prior  : num [1:2] 1 1
#>   .. ..$ nchains: num 2
#>   .. ..$ burnin : num 10000
#>   .. ..$ update : num 10000
#>   .. ..$ inits  : NULL
#>   ..@ model      :Class 'prevModel'  chr [1:7] "model {" "x ~ dbin(AP, n)" "AP <- SE * TP + (1 - SP) * (1 - TP)" "SE ~ dunif(0.6, 1)" ...
#>   ..@ mcmc       :List of 3
#>   .. ..$ TP:List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.125 0.1236 0.1076 0.0944 0.133 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.111 0.0971 0.1033 0.1291 0.1516 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   .. ..$ SE:List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.673 0.696 0.788 0.693 0.629 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.825 0.834 0.88 0.806 0.835 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   .. ..$ SP:List of 2
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.877 0.87 0.868 0.863 0.861 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..$ :Class 'mcmc'  atomic [1:10000] 0.884 0.875 0.917 0.89 0.927 ...
#>   .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   .. .. ..- attr(*, "class")= chr "mcmc.list"
#>   ..@ diagnostics:List of 2
#>   .. ..$ DIC:List of 3
#>   .. .. ..$ deviance:Class 'mcarray'  Named num 7.58
#>   .. .. .. .. ..- attr(*, "names")= chr "x"
#>   .. .. ..$ penalty :Class 'mcarray'  Named num 0.968
#>   .. .. .. .. ..- 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] "SE" "SP" "TP"
#>   .. .. .. .. ..$ : chr [1:2] "Point est." "Upper C.I."
#>   .. .. ..$ mpsrf: num 1
#>   .. .. ..- attr(*, "class")= chr "gelman.diag"
## Every slot can be accessed using the '@' operator
slotNames(TP)
#> [1] "par"         "model"       "mcmc"        "diagnostics"
## .. input parameters
TP@par
#> $x
#> [1] 142
#> 
#> $n
#> [1] 742
#> 
#> $SE
#> $SE$d
#> [1] "uniform"
#> 
#> $SE$p
#> [1] 0.6 1.0
#> 
#> 
#> $SP
#> $SP$d
#> [1] "uniform"
#> 
#> $SP$p
#> [1] 0.75 1.00
#> 
#> 
#> $prior
#> [1] 1 1
#> 
#> $nchains
#> [1] 2
#> 
#> $burnin
#> [1] 10000
#> 
#> $update
#> [1] 10000
#> 
#> $inits
#> NULL
## .. fitted model
TP@model
#> model {
#>   x ~ dbin(AP, n)
#>   AP <- SE * TP + (1 - SP) * (1 - TP)
#>   SE ~ dunif(0.6, 1)
#>   SP ~ dunif(0.75, 1)
#>   TP ~ dbeta(1, 1)
#> }
## .. simulated TP, SE, SP
str(TP@mcmc)
#> List of 3
#>  $ TP:List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.125 0.1236 0.1076 0.0944 0.133 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.111 0.0971 0.1033 0.1291 0.1516 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
#>  $ SE:List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.673 0.696 0.788 0.693 0.629 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.825 0.834 0.88 0.806 0.835 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
#>  $ SP:List of 2
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.877 0.87 0.868 0.863 0.861 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..$ :Class 'mcmc'  atomic [1:10000] 0.884 0.875 0.917 0.89 0.927 ...
#>   .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#>   ..- attr(*, "class")= chr "mcmc.list"
## .. convert simulated TP, SE, SP to matrix
head(as.matrix(TP))
#>              TP        SE        SP
#> [1,] 0.12497886 0.6729475 0.8774353
#> [2,] 0.12357302 0.6959699 0.8699901
#> [3,] 0.10761131 0.7880032 0.8683324
#> [4,] 0.09438967 0.6927515 0.8632794
#> [5,] 0.13304616 0.6286654 0.8607214
#> [6,] 0.12531370 0.7360327 0.8846877
## .. DIC and BGR (and bayesP)
TP@diagnostics
#> $DIC
#> Mean deviance:  7.576 
#> penalty 0.9684 
#> Penalized deviance: 8.544 
#> 
#> $BGR
#> Potential scale reduction factors:
#> 
#>    Point est. Upper C.I.
#> SE          1          1
#> SP          1          1
#> TP          1          1
#> 
#> Multivariate psrf
#> 
#> 1
## TP@mcmc elements inherit from class 'mcmc.list' in coda package
## List all available methods for this class
methods(class = "mcmc.list")
#>  [1] [             acfplot       as.array      as.matrix     as.mcmc      
#>  [6] autocorr.diag batchSE       coerce        densityplot   end          
#> [11] head          HPDinterval   initialize    plot          qqmath       
#> [16] rejectionRate show          slotsFromS3   start         summary      
#> [21] tail          thin          time          window        xyplot       
#> see '?methods' for accessing help and source code