View source: R/truePrevMulti-main.R
truePrevMulti2 | 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 ≥ 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 |
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.