the prevalence package
tools for prevalence assessment studies.
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 method
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