truePrevMulti2 {prevalence} | R Documentation |
Bayesian estimation of true prevalence from apparent prevalence obtained by applying multiple tests to individual samples. truePrevMulti2
implements and extends the approach described by Dendukuri and Joseph (2001), which uses a multinomial distribution to model observed test results,
and in which conditional dependence between tests is modelled through
covariances.
truePrevMulti2(x, n, prior, nchains = 2, burnin = 10000, update = 10000,
verbose = FALSE)
x |
Vector of apparent test results; see 'Details' below |
n |
The total sample size |
prior |
The prior distributions; see 'Details' below |
nchains |
The number of chains used in the estimation process; must be |
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 |
truePrevMulti2
calls on JAGS via the rjags package to estimate true prevalence from apparent prevalence in a Bayesian framework. truePrevMulti2
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
-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
.
Argument prior
consists of prior distributions for:
True Prevalence: TP
SEnsitivity of each individual test: vector SE
SPecificity of each individual test: vector SP
Conditional covariance of all possible test combinations given a truly positive disease status: vector a
Conditional covariance of all possible test combinations given a truly negative disease status: vector b
To see how prior
is defined for the number of tests h
at hand, use define_prior2
.
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 prior
should be given between curly brackets (i.e., {}
), separated by line breaks. Priors 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, priors 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.
An object of class prev
.
Markov chain Monte Carlo sampling in truePrevMulti2
is performed by JAGS (Just Another Gibbs Sampler) through the rjags package. JAGS can be downloaded from https://mcmc-jags.sourceforge.io/.
Brecht Devleesschauwer <brechtdv@gmail.com>
Dendukuri N, Joseph L (2001) Bayesian approaches to modeling the conditional dependence between multiple diagnostic tests. Biometrics 57:158-167
define_x
: how to define the vector of apparent test results x
define_prior2
: how to define prior
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
truePrev
: estimate true prevalence from apparent prevalence obtained by testing individual samples with a single test
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
## Not run:
## ===================================================== ##
## 2-TEST EXAMPLE: Strongyloides ##
## ----------------------------------------------------- ##
## Two tests were performed on 162 humans ##
## -> T1 = stool examination ##
## -> T2 = serology test ##
## Expert opinion generated the following priors: ##
## -> SE1 ~ dbeta( 4.44, 13.31) ##
## -> SP1 ~ dbeta(71.25, 3.75) ##
## -> SE2 ~ dbeta(21.96, 5.49) ##
## -> SP2 ~ dbeta( 4.10, 1.76) ##
## The following results were obtained: ##
## -> 38 samples T1+,T2+ ##
## -> 2 samples T1+,T2- ##
## -> 87 samples T1-,T2+ ##
## -> 35 samples T1-,T2- ##
## ===================================================== ##
## how is the 2-test model defined?
define_x(2)
define_prior2(2)
## fit Strongyloides 2-test model
## a first model assumes conditional independence
## -> set covariance terms to zero
strongy_indep <-
truePrevMulti2(
x = c(38, 2, 87, 35),
n = 162,
prior = {
TP ~ dbeta(1, 1)
SE[1] ~ dbeta( 4.44, 13.31)
SP[1] ~ dbeta(71.25, 3.75)
SE[2] ~ dbeta(21.96, 5.49)
SP[2] ~ dbeta( 4.10, 1.76)
a[1] <- 0
b[1] <- 0
})
## show model results
strongy_indep
## fit same model using 'list-style'
strongy_indep <-
truePrevMulti2(
x = c(38, 2, 87, 35),
n = 162,
prior =
list(
TP = list(dist = "beta", alpha = 1, beta = 1),
SE1 = list(dist = "beta", alpha = 4.44, beta = 13.31),
SP1 = list(dist = "beta", alpha = 71.25, beta = 3.75),
SE2 = list(dist = "beta", alpha = 21.96, beta = 5.49),
SP2 = list(dist = "beta", alpha = 4.10, beta = 1.76),
a1 = 0,
b1 = 0
)
)
## show model results
strongy_indep
## fit Strongyloides 2-test model
## a second model allows for conditional dependence
## -> a[1] is the covariance between T1 and T2, given D+
## -> b[1] is the covariance between T1 and T2, given D-
## -> a[1] and b[1] can range between +/- 2^-h, ie, (-.25, .25)
strongy <-
truePrevMulti2(
x = c(38, 2, 87, 35),
n = 162,
prior = {
TP ~ dbeta(1, 1)
SE[1] ~ dbeta( 4.44, 13.31)
SP[1] ~ dbeta(71.25, 3.75)
SE[2] ~ dbeta(21.96, 5.49)
SP[2] ~ dbeta( 4.10, 1.76)
a[1] ~ dunif(-0.25, 0.25)
b[1] ~ dunif(-0.25, 0.25)
})
## explore model structure
str(strongy) # overall structure
str(strongy@par) # structure of slot 'par'
str(strongy@mcmc) # structure of slot 'mcmc'
strongy@model # fitted model
strongy@diagnostics # DIC, BGR and Bayes-P values
## standard methods
print(strongy)
summary(strongy)
par(mfrow = c(2, 2))
plot(strongy) # shows plots of TP by default
plot(strongy, "SE[1]") # same plots for SE1
plot(strongy, "SE[2]") # same plots for SE2
plot(strongy, "SP[1]") # same plots for SP1
plot(strongy, "SP[2]") # same plots for SP2
plot(strongy, "a[1]") # same plots for a[1]
plot(strongy, "b[1]") # same plots for b[1]
## coda plots of all parameters
par(mfrow = c(2, 4)); densplot(strongy, col = "red")
par(mfrow = c(2, 4)); traceplot(strongy)
par(mfrow = c(2, 4)); gelman.plot(strongy)
par(mfrow = c(2, 4)); autocorr.plot(strongy)
## End(Not run)