The ptb R package is a work in progress to implement functions, under a Bayesian framework, for prevalence estimation, diagnostic test evaluation and prior elicitation. The ptb runs with JAGS (http://sourceforge.net/projects/mcmc-jags/) as backend. The functions implemented so far are generalizations of code presented by professor Ian Garner in a course he offered in São Paulo, Brazil. As ptb functions return R objects, other nice packages can be used for model diagnostics.
After installing JAGS (http://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/), use the following commands:
library(devtools) install_github("leb-fmvz-usp/ptb")
knitr::opts_chunk$set(echo = TRUE)
library(ptb); library(coda); library(ggmcmc)
The function ElicitBeta
can be used to calculate the parameters a and b of a Beta distribution. The function takes as imputs the elicited mode for the variable of interest and the maximum or minimum elicited value for that variable. A confidence about this maximum or minimum must be provided too. The function PlotElicitedPrior
plots the resulting ditribution.
(se_prior <- ElicitBeta(mode = 0.3, maximum = 0.5, confidence = 0.95)) PlotElicitedPrior(se_prior)
The argument summary
provides summary statistics and the argument quantiles
allows the specification of specific quantiles.
ElicitBeta(mode = 0.3, maximum = 0.5, confidence = 0.95, summary = TRUE) ElicitBeta(mode = 0.3, maximum = 0.5, confidence = 0.95, summary = TRUE, quantiles = c(0.1, 0.9))
To compile and run the models, we always need to define the data, the priors and the parameters to be monitored.
By default, three parallel chains are run, the “burn in” is equal to 1e3, starting values are automatically generated, the number of effective iterations is equal to 1e4 and chains are not “thinned”. For details, see the help pages.
# help(OneTestOnePopBM) # Data dataset <- list(pop_size = 91, positives = 1) # Initial conditions for chains (optional) inits <- list(list(true_prev = 0.05, se = 0.8, sp = 0.9), list(true_prev = 0.02, se = 0.3, sp = 0.7), list(true_prev = 0.09, se = 0.1, sp = 0.5)) # Priors priors <- c(true_prev_a = 1, true_prev_b = 1, se_a = 6.28, se_b = 13.32, sp_a = 212.12, sp_b = 3.13) # Prevalence estimate prev_est <- OneTestOnePopBM(dataset = dataset, inits = inits, priors = priors, pars = c('true_prev', 'ppv', 'npv')) summary(prev_est)
# help(OneTestOnePopBMM) # Data dataset <- list(pop_size = 91, positives = 1) # Priors priors <- list(true_prev_wph_a = 1.8, true_prev_wph_b = 26.74, se_a = 6.28, se_b = 13.32, sp_a = 212.12, sp_b = 3.13, prev_h = 0.1) # Prevalence estimates prev_est <- OneTestOnePopBMM(dataset = dataset, priors = priors, n_iter = 3e3, pars = c('true_prev', 'true_prev_wph', 'prev_h')) summary(prev_est)
# help(DiffBetweenEstimates) # Priors priors <- c(true_prev_a = 1, true_prev_b = 1, se_a = 6.28, se_b = 13.32, sp_a = 212.12, sp_b = 3.13) # First estimate dataset1 <- list(pop_size = 100, positives = 5) prev_est1 <- OneTestOnePopBM(dataset = as.list(dataset1), n_iter = 3e3, priors = priors, pars = "true_prev", burn_in = 5e2) # Second estimate dataset2 <- list(pop_size = 91, positives = 1) prev_est2 <- OneTestOnePopBM(dataset = as.list(dataset2), n_iter = 3e3, priors = priors, pars = "true_prev", burn_in = 5e2) # Estimated difference. diffs <- DiffBetweenEstimates(list(prev_est1, prev_est2)) summary(diffs)
The model outputs are ready for diagnostics.
gelman.diag(prev_est) gelman.plot(prev_est) gg_res <- ggs(prev_est) ggs_traceplot(gg_res) ggs_density(gg_res) ggs_histogram(gg_res, bins = 100) ggs_compare_partial(gg_res) ggs_running(gg_res) ggs_autocorrelation(gg_res)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.