In leb-fmvz-usp/ptb: Prevalence Estimation and Diagnostic Test Evaluation in a Bayesian Framework

ptb: an R package for prevalence estimation and diagnostic test evaluation in a Bayesian framework

Oswaldo Santos Baquero

contact: [e-mail]

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.

Installation

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)

Load packages

library(ptb); library(coda); library(ggmcmc)

Fit a Beta distribution from elicited information

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))

Models

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.

One test and one population binomial model

# 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)

One test and one population binomial mixture model

# 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)

Difference between estimates

# 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)

Diagnostics

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)

leb-fmvz-usp/ptb documentation built on May 30, 2019, 3:43 p.m.