the prevalence package
tools for prevalence assessment studies.
tools for prevalence assessment studies.
Estimate true prevalence from individuals samples using multiple tests – covariance scheme
Description
Bayesian estimation of true prevalence from apparent prevalence obtained by applying multiple tests to individual samples. truePrevMulti2 implements 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.
Usage
truePrevMulti2(x, n, prior, nchains = 2, burnin = 10000, update = 10000,
verbose = FALSE)Arguments
x
|
Vector of apparent test results; see details |
n
|
The total sample size |
prior
|
The prior distribution for theta; see details
|
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
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-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.
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 n, 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 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.
Alternatively, priors can be specified in a named list().
More info on specifying distributions is available here.
Value
A prev-class object.
Note
Markov chain Monte Carlo sampling in truePrevMulti2 is performed by JAGS (Just Another Gibbs Sampler) through the rjags package. JAGS can be downloaded from http://sourceforge.net/projects/mcmc-jags/.
References
- Dendukuri N, Joseph L (2001) Bayesian approaches to modeling the conditional dependence between multiple diagnostic tests. Biometrics 57:158-167. http://dx.doi.org/10.1111/j.0006-341X.2001.00158.x
See Also
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 with a single test
betaPERT: calculate the parameters of a Beta-PERT distribution
betaExpert: calculate the parameters of a Beta distribution based on expert opinion
Examples
## ===================================================== ##
## 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)
#> 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_prior2(2)
#> Covariance scheme
#> Definition of the prior for 2 tests:
#> TP : True Prevalence
#> SE[1] : Sensitity T1
#> SE[2] : Sensitity T2
#> SP[1] : Specificity T1
#> SP[2] : Specificity T2
#> a[1] : Covariance(T1,T2|D+)
#> b[1] : Covariance(T1,T2|D-)
## 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
#> mean median mode sd 2.5% 97.5%
#> TP 0.736 0.752 0.774 0.086 0.526 0.865
#> SE[1] 0.324 0.317 0.313 0.047 0.251 0.437
#> SE[2] 0.914 0.917 0.927 0.030 0.845 0.964
#> SP[1] 0.963 0.966 0.971 0.018 0.921 0.990
#> SP[2] 0.703 0.712 0.744 0.152 0.405 0.957
#> a[1] 0.000 0.000 0.000 0.000 0.000 0.000
#> b[1] 0.000 0.000 0.000 0.000 0.000 0.000
#>
#> Multivariate BGR statistic = 1.0005
#> BGR values substantially above 1 indicate lack of convergence
## 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
#> mean median mode sd 2.5% 97.5%
#> TP 0.740 0.755 0.777 0.084 0.533 0.864
#> SE[1] 0.322 0.316 0.310 0.045 0.252 0.427
#> SE[2] 0.914 0.917 0.916 0.030 0.847 0.964
#> SP[1] 0.963 0.966 0.972 0.018 0.921 0.990
#> SP[2] 0.710 0.720 0.745 0.150 0.408 0.959
#> a[1] 0.000 0.000 0.000 0.000 0.000 0.000
#> b[1] 0.000 0.000 0.000 0.000 0.000 0.000
#>
#> Multivariate BGR statistic = 1.0024
#> BGR values substantially above 1 indicate lack of convergence
## 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
#> Formal class 'prev' [package "prevalence"] with 4 slots
#> ..@ par :List of 7
#> .. ..$ x : num [1:4] 38 2 87 35
#> .. ..$ n : num 162
#> .. ..$ prior :List of 7
#> .. .. ..$ :List of 2
#> .. .. .. ..$ : chr "TP"
#> .. .. .. ..$ :List of 2
#> .. .. .. .. ..$ d: chr "beta"
#> .. .. .. .. ..$ p: num [1:2] 1 1
#> .. .. ..$ :List of 2
#> .. .. .. ..$ : chr "SE[1]"
#> .. .. .. ..$ :List of 2
#> .. .. .. .. ..$ d: chr "beta"
#> .. .. .. .. ..$ p: num [1:2] 4.44 13.31
#> .. .. ..$ :List of 2
#> .. .. .. ..$ : chr "SE[2]"
#> .. .. .. ..$ :List of 2
#> .. .. .. .. ..$ d: chr "beta"
#> .. .. .. .. ..$ p: num [1:2] 21.96 5.49
#> .. .. ..$ :List of 2
#> .. .. .. ..$ : chr "SP[1]"
#> .. .. .. ..$ :List of 2
#> .. .. .. .. ..$ d: chr "beta"
#> .. .. .. .. ..$ p: num [1:2] 71.25 3.75
#> .. .. ..$ :List of 2
#> .. .. .. ..$ : chr "SP[2]"
#> .. .. .. ..$ :List of 2
#> .. .. .. .. ..$ d: chr "beta"
#> .. .. .. .. ..$ p: num [1:2] 4.1 1.76
#> .. .. ..$ :List of 2
#> .. .. .. ..$ : chr "a[1]"
#> .. .. .. ..$ :List of 2
#> .. .. .. .. ..$ d: chr "uniform"
#> .. .. .. .. ..$ p: num [1:2] -0.25 0.25
#> .. .. ..$ :List of 2
#> .. .. .. ..$ : chr "b[1]"
#> .. .. .. ..$ :List of 2
#> .. .. .. .. ..$ d: chr "uniform"
#> .. .. .. .. ..$ p: num [1:2] -0.25 0.25
#> .. ..$ nchains: num 2
#> .. ..$ burnin : num 10000
#> .. ..$ update : num 10000
#> .. ..$ inits : NULL
#> ..@ model :Class 'prevModel' chr [1:47] "model {" "x[1:4] ~ dmulti(AP[1:4], n)" "" "prob_se[1] <- SE[1] * SE[2] + a[1]" ...
#> ..@ mcmc :List of 8
#> .. ..$ TP :List of 2
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.606 0.577 0.573 0.494 0.638 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.791 0.832 0.862 0.813 0.829 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..- attr(*, "class")= chr "mcmc.list"
#> .. ..$ SE[1] :List of 2
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.387 0.37 0.354 0.383 0.378 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.263 0.268 0.29 0.303 0.271 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..- attr(*, "class")= chr "mcmc.list"
#> .. ..$ SE[2] :List of 2
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.834 0.836 0.848 0.874 0.856 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.864 0.866 0.872 0.87 0.843 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..- attr(*, "class")= chr "mcmc.list"
#> .. ..$ SP[1] :List of 2
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.908 0.912 0.916 0.932 0.925 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.906 0.897 0.912 0.936 0.921 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..- attr(*, "class")= chr "mcmc.list"
#> .. ..$ SP[2] :List of 2
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.431 0.421 0.413 0.401 0.44 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.75 0.655 0.645 0.601 0.712 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..- attr(*, "class")= chr "mcmc.list"
#> .. ..$ a[1] :List of 2
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.0403 0.0367 0.0301 0.0307 0.0315 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.0352 0.0318 0.0278 0.0307 0.0337 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..- attr(*, "class")= chr "mcmc.list"
#> .. ..$ b[1] :List of 2
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.02334 0.03267 0.01449 0.00939 0.01021 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..$ :Class 'mcmc' atomic [1:10000] 0.0044 -0.0304 -0.0144 -0.011 0.0223 ...
#> .. .. .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> .. .. ..- attr(*, "class")= chr "mcmc.list"
#> .. ..$ bayesP: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"
#> ..@ diagnostics:List of 3
#> .. ..$ DIC :List of 3
#> .. .. ..$ deviance:Class 'mcarray' Named num [1:26] 14.4 14.4 0 0 0 ...
#> .. .. .. .. ..- attr(*, "names")= chr [1:26] "x" "x2" "O6[1]" "O5[1]" ...
#> .. .. ..$ penalty :Class 'mcarray' Named num [1:26] 64749 64749 0 0 0 ...
#> .. .. .. .. ..- attr(*, "names")= chr [1:26] "x" "x2" "O6[1]" "O5[1]" ...
#> .. .. ..$ type : chr "pD"
#> .. .. ..- attr(*, "class")= chr "dic"
#> .. ..$ BGR :List of 2
#> .. .. ..$ psrf : num [1:7, 1:2] 1 1 1 1 1 ...
#> .. .. .. ..- attr(*, "dimnames")=List of 2
#> .. .. .. .. ..$ : chr [1:7] "SE[1]" "SE[2]" "SP[1]" "SP[2]" ...
#> .. .. .. .. ..$ : chr [1:2] "Point est." "Upper C.I."
#> .. .. ..$ mpsrf: num 1
#> .. .. ..- attr(*, "class")= chr "gelman.diag"
#> .. ..$ bayesP: num 1
str(strongy@par) # structure of slot 'par'
#> List of 7
#> $ x : num [1:4] 38 2 87 35
#> $ n : num 162
#> $ prior :List of 7
#> ..$ :List of 2
#> .. ..$ : chr "TP"
#> .. ..$ :List of 2
#> .. .. ..$ d: chr "beta"
#> .. .. ..$ p: num [1:2] 1 1
#> ..$ :List of 2
#> .. ..$ : chr "SE[1]"
#> .. ..$ :List of 2
#> .. .. ..$ d: chr "beta"
#> .. .. ..$ p: num [1:2] 4.44 13.31
#> ..$ :List of 2
#> .. ..$ : chr "SE[2]"
#> .. ..$ :List of 2
#> .. .. ..$ d: chr "beta"
#> .. .. ..$ p: num [1:2] 21.96 5.49
#> ..$ :List of 2
#> .. ..$ : chr "SP[1]"
#> .. ..$ :List of 2
#> .. .. ..$ d: chr "beta"
#> .. .. ..$ p: num [1:2] 71.25 3.75
#> ..$ :List of 2
#> .. ..$ : chr "SP[2]"
#> .. ..$ :List of 2
#> .. .. ..$ d: chr "beta"
#> .. .. ..$ p: num [1:2] 4.1 1.76
#> ..$ :List of 2
#> .. ..$ : chr "a[1]"
#> .. ..$ :List of 2
#> .. .. ..$ d: chr "uniform"
#> .. .. ..$ p: num [1:2] -0.25 0.25
#> ..$ :List of 2
#> .. ..$ : chr "b[1]"
#> .. ..$ :List of 2
#> .. .. ..$ d: chr "uniform"
#> .. .. ..$ p: num [1:2] -0.25 0.25
#> $ nchains: num 2
#> $ burnin : num 10000
#> $ update : num 10000
#> $ inits : NULL
str(strongy@mcmc) # structure of slot 'mcmc'
#> List of 8
#> $ TP :List of 2
#> ..$ :Class 'mcmc' atomic [1:10000] 0.606 0.577 0.573 0.494 0.638 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..$ :Class 'mcmc' atomic [1:10000] 0.791 0.832 0.862 0.813 0.829 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ SE[1] :List of 2
#> ..$ :Class 'mcmc' atomic [1:10000] 0.387 0.37 0.354 0.383 0.378 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..$ :Class 'mcmc' atomic [1:10000] 0.263 0.268 0.29 0.303 0.271 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ SE[2] :List of 2
#> ..$ :Class 'mcmc' atomic [1:10000] 0.834 0.836 0.848 0.874 0.856 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..$ :Class 'mcmc' atomic [1:10000] 0.864 0.866 0.872 0.87 0.843 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ SP[1] :List of 2
#> ..$ :Class 'mcmc' atomic [1:10000] 0.908 0.912 0.916 0.932 0.925 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..$ :Class 'mcmc' atomic [1:10000] 0.906 0.897 0.912 0.936 0.921 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ SP[2] :List of 2
#> ..$ :Class 'mcmc' atomic [1:10000] 0.431 0.421 0.413 0.401 0.44 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..$ :Class 'mcmc' atomic [1:10000] 0.75 0.655 0.645 0.601 0.712 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ a[1] :List of 2
#> ..$ :Class 'mcmc' atomic [1:10000] 0.0403 0.0367 0.0301 0.0307 0.0315 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..$ :Class 'mcmc' atomic [1:10000] 0.0352 0.0318 0.0278 0.0307 0.0337 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ b[1] :List of 2
#> ..$ :Class 'mcmc' atomic [1:10000] 0.02334 0.03267 0.01449 0.00939 0.01021 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..$ :Class 'mcmc' atomic [1:10000] 0.0044 -0.0304 -0.0144 -0.011 0.0223 ...
#> .. .. ..- attr(*, "mcpar")= num [1:3] 11001 21000 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ bayesP: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"
strongy@model # fitted model
#> model {
#> x[1:4] ~ dmulti(AP[1:4], n)
#>
#> prob_se[1] <- SE[1] * SE[2] + a[1]
#> prob_se[2] <- SE[1] * (1 - SE[2]) - a[1]
#> prob_se[3] <- (1 - SE[1]) * SE[2] - a[1]
#> prob_se[4] <- (1 - SE[1]) * (1 - SE[2]) + a[1]
#>
#> prob_sp[4] <- SP[1] * SP[2] + b[1]
#> prob_sp[3] <- SP[1] * (1 - SP[2]) - b[1]
#> prob_sp[2] <- (1 - SP[1]) * SP[2] - b[1]
#> prob_sp[1] <- (1 - SP[1]) * (1 - SP[2]) + b[1]
#>
#> for (i in 1:4) {
#> AP[i] <- TP * prob_se[i] + (1 - TP) * prob_sp[i]
#>
#> constraint1[i] <- step(AP[i])
#> O1[i] ~ dbern(constraint1[i])
#> constraint2[i] <- step(AP[i] - 1)
#> O2[i] ~ dbern(constraint2[i])
#> constraint3[i] <- step(prob_se[i])
#> O3[i] ~ dbern(constraint3[i])
#> constraint4[i] <- step(prob_se[i] - 1)
#> O4[i] ~ dbern(constraint4[i])
#> constraint5[i] <- step(prob_sp[i])
#> O5[i] ~ dbern(constraint5[i])
#> constraint6[i] <- step(prob_sp[i] - 1)
#> O6[i] ~ dbern(constraint6[i])
#> }
#>
#> TP ~ dbeta(1, 1)
#> SE[1] ~ dbeta(4.44, 13.31)
#> SE[2] ~ dbeta(21.96, 5.49)
#> SP[1] ~ dbeta(71.25, 3.75)
#> SP[2] ~ dbeta(4.1, 1.76)
#> a[1] ~ dunif(-0.25, 0.25)
#> b[1] ~ dunif(-0.25, 0.25)
#>
#> 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 <- sum(d1[])
#> Gt <- sum(d2[])
#> bayesP <- step(G0 - Gt)
#> }
strongy@diagnostics # DIC, BGR and Bayes-P values
#> $DIC
#> Mean deviance: 28.75
#> penalty 129498
#> Penalized deviance: 129527
#>
#> $BGR
#> Potential scale reduction factors:
#>
#> Point est. Upper C.I.
#> SE[1] 1 1.01
#> SE[2] 1 1.00
#> SP[1] 1 1.00
#> SP[2] 1 1.01
#> TP 1 1.00
#> a 1 1.01
#> b 1 1.01
#>
#> Multivariate psrf
#>
#> 1
#>
#> $bayesP
#> [1] 1
## standard methods
print(strongy)
#> mean median mode sd 2.5% 97.5%
#> TP 0.842 0.861 0.889 0.109 0.577 0.991
#> SE[1] 0.284 0.278 0.270 0.042 0.219 0.387
#> SE[2] 0.831 0.828 0.821 0.046 0.750 0.923
#> SP[1] 0.938 0.942 0.949 0.027 0.876 0.980
#> SP[2] 0.636 0.648 0.666 0.186 0.278 0.937
#> a[1] 0.033 0.034 0.037 0.014 0.005 0.059
#> b[1] 0.011 0.009 0.003 0.022 -0.028 0.061
#>
#> Multivariate BGR statistic = 1.0028
#> BGR values substantially above 1 indicate lack of convergence
summary(strongy)
#> $TP
#> mean median mode sd var 2.5%
#> chain 1 0.8436507 0.8610896 0.8909832 0.1056432 0.01116048 0.5935869
#> chain 2 0.8399529 0.8601476 0.8884333 0.1118842 0.01251807 0.5600494
#> all chains 0.8418018 0.8606750 0.8892378 0.1088214 0.01184210 0.5771182
#> 97.5% samples
#> chain 1 0.9910774 10000
#> chain 2 0.9907360 10000
#> all chains 0.9909564 20000
#>
#> $`SE[1]`
#> mean median mode sd var 2.5%
#> chain 1 0.2832589 0.2778963 0.2678959 0.04053746 0.001643286 0.2185224
#> chain 2 0.2856656 0.2788679 0.2719301 0.04395325 0.001931888 0.2201134
#> all chains 0.2844622 0.2784071 0.2695707 0.04229593 0.001788946 0.2190952
#> 97.5% samples
#> chain 1 0.3793029 10000
#> chain 2 0.3942700 10000
#> all chains 0.3871296 20000
#>
#> $`SE[2]`
#> mean median mode sd var 2.5%
#> chain 1 0.8318448 0.8290919 0.8248206 0.04586113 0.002103243 0.7505227
#> chain 2 0.8298677 0.8267019 0.8155388 0.04566812 0.002085577 0.7488405
#> all chains 0.8308563 0.8279175 0.8211626 0.04577426 0.002095283 0.7496565
#> 97.5% samples
#> chain 1 0.9242968 10000
#> chain 2 0.9207437 10000
#> all chains 0.9228179 20000
#>
#> $`SP[1]`
#> mean median mode sd var 2.5%
#> chain 1 0.9380062 0.9418135 0.9493370 0.02654011 0.0007043775 0.8769372
#> chain 2 0.9386962 0.9428199 0.9500813 0.02744737 0.0007533582 0.8745201
#> all chains 0.9383512 0.9423324 0.9493613 0.02699908 0.0007289504 0.8758966
#> 97.5% samples
#> chain 1 0.9787722 10000
#> chain 2 0.9802919 10000
#> all chains 0.9796821 20000
#>
#> $`SP[2]`
#> mean median mode sd var 2.5%
#> chain 1 0.6411893 0.6531746 0.6486065 0.1830358 0.03350210 0.2810302
#> chain 2 0.6303553 0.6442841 0.7094148 0.1880421 0.03535984 0.2734552
#> all chains 0.6357723 0.6484476 0.6660765 0.1856303 0.03445859 0.2776680
#> 97.5% samples
#> chain 1 0.9348416 10000
#> chain 2 0.9384042 10000
#> all chains 0.9365680 20000
#>
#> $`a[1]`
#> mean median mode sd var
#> chain 1 0.03281355 0.03349853 0.03569904 0.01354784 0.0001835440
#> chain 2 0.03367770 0.03415930 0.03698381 0.01407399 0.0001980773
#> all chains 0.03324563 0.03383647 0.03707119 0.01381983 0.0001909878
#> 2.5% 97.5% samples
#> chain 1 0.004390537 0.05811966 10000
#> chain 2 0.004786837 0.06002443 10000
#> all chains 0.004588747 0.05927007 20000
#>
#> $`b[1]`
#> mean median mode sd var
#> chain 1 0.01183860 0.009658547 0.004278029 0.02248901 0.0005057554
#> chain 2 0.01040055 0.008599015 0.002992741 0.02175976 0.0004734874
#> all chains 0.01111957 0.009129480 0.002841077 0.02213852 0.0004901139
#> 2.5% 97.5% samples
#> chain 1 -0.02838421 0.06168614 10000
#> chain 2 -0.02858627 0.05925678 10000
#> all chains -0.02839380 0.06059105 20000
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)
