bamdit (version 3.4.4)

metadiag: Bayesian Meta-Analysis of diagnostic test data

Description

This function performers a Bayesian meta-analysis of diagnostic test data by fitting a bivariate random effects model. The number of true positives and false positives are modeled with two conditional Binomial distributions and the random-effects are based on a bivariate scale mixture of Normals. Computations are done by calling JAGS (Just Another Gibbs Sampler) to perform MCMC (Markov Chain Monte Carlo) sampling and returning an object of the class mcmc.list.

Usage

metadiag(
  data,
  two.by.two = FALSE,
  re = "normal",
  re.model = "DS",
  link = "logit",
  mean.mu.D = 0,
  mean.mu.S = 0,
  sd.mu.D = 1,
  sd.mu.S = 1,
  sigma.D.upper = 10,
  sigma.S.upper = 10,
  mean.Fisher.rho = 0,
  sd.Fisher.rho = 1/sqrt(2),
  df = 4,
  df.estimate = FALSE,
  df.lower = 3,
  df.upper = 20,
  split.w = FALSE,
  n.1.new = 50,
  n.2.new = 50,
  nr.chains = 2,
  nr.iterations = 10000,
  nr.adapt = 1000,
  nr.burnin = 1000,
  nr.thin = 1,
  be.quiet = FALSE,
  r2jags = TRUE
)

Value

This function returns an object of the class metadiag. This object contains the MCMC output of each parameter and hyper-parameter in the model, the data frame used for fitting the model, the link function, type of random effects distribution and the splitting information for conflict of evidence analysis.

The results of the object of the class metadiag can be extracted with R2jags or with rjags. In addition a summary, a print and a plot functions are implemented for this type of object.

Arguments

data

Either a data frame with at least 4 columns containing the true positives (tp), number of patients with disease (n1), false positives (fp), number of patients without disease (n2), or for two.by.two = TRUE a data frame where each line contains the diagnostic results as a two by two table, where the column names are: TP, FP, TN, FN.

two.by.two

If TRUE indicates that the diagnostic results are given as: TP, FP, TN, FN.

re

Random effects distribution for the resulting model. Possible values are normal for bivariate random effects and sm for scale mixtures.

re.model

If re.model = "DS" indicates that the sum and differences of TPR and FPR are modeled as random effects and re.model = "SeSp" indicates that the Sensitivity and Specificity are modeled as ranodm effects. The default value is re.model = "DS".

link

The link function used in the model. Possible values are logit, cloglog probit.

mean.mu.D

prior Mean of D, default value is 0.

mean.mu.S

prior Mean of S, default value is 0.

sd.mu.D

prior Standard deviation of D, default value is 1 (the prior of mu.D is a logistic distribution).

sd.mu.S

prior Standard deviation of S, default value is 1 (the prior of mu.S is a logistic distribution).

sigma.D.upper

Upper bound of the uniform prior of sigma.S, default value is 10.

sigma.S.upper

Upper bound of the uniform prior of sigma.S, default value is 10.

mean.Fisher.rho

Mean of rho in the Fisher scale default value is 0.

sd.Fisher.rho

Standard deviation of rho in the Fisher scale, default value is 1/sqrt(2).

df

If de.estimate = FALSE, then df is the degrees of freedom for the scale mixture distribution, default value is 4.

df.estimate

Estimate the posterior of df. The defualt value is FALSE.

df.lower

Lower bound of the prior of df. The defulat value is 3.

df.upper

Upper bound of the prior of df. The defulat value is 30.

split.w

Split the w parameter in two independent weights one for each random effect. The default value is FALSE.

n.1.new

Number of patients with disease in a predictive study default is 50.

n.2.new

Number of patients with non-disease in a predictive study default is 50.

nr.chains

Number of chains for the MCMC computations, default 5.

nr.iterations

Number of iterations after adapting the MCMC, default is 10000. Some models may need more iterations.

nr.adapt

Number of iterations in the adaptation process, default is 1000. Some models may need more iterations during adaptation.

nr.burnin

Number of iteration discard for burn-in period, default is 1000. Some models may need a longer burnin period.

nr.thin

Thinning rate, it must be a positive integer, the default value 1.

be.quiet

Do not print warning message if the model does not adapt default value is FALSE. If you are not sure about the adaptation period choose be.quiet=TRUE.

r2jags

Which interface is used to link R to JAGS (rjags and R2jags) default value is R2Jags TRUE.

Details

Installation of JAGS: It is important to note that R 3.3.0 introduced a major change in the use of toolchain for Windows. This new toolchain is incompatible with older packages written in C++. As a consequence, if the installed version of JAGS does not match the R installation, then the rjags package will spontaneously crash. Therefore, if a user works with R version >= 3.3.0, then JAGS must be installed with the installation program JAGS-4.2.0-Rtools33.exe. For users who continue using R 3.2.4 or an earlier version, the installation program for JAGS is the default installer JAGS-4.2.0.exe.

References

Verde P. E. (2010). Meta-analysis of diagnostic test data: A bivariate Bayesian modeling approach. Statistics in Medicine. 29(30):3088-102. doi: 10.1002/sim.4055.

Verde P. E. (2018). bamdit: An R Package for Bayesian Meta-Analysis of Diagnostic Test Data. Journal of Statisticsl Software. Volume 86, issue 10, pages 1--32.

Examples

if (FALSE) {

# Example: data from Glas et al. (2003).....................................
library(bamdit)
data("glas")
glas.t <- glas[glas$marker == "Telomerase", 1:4]

glas.t <- glas[glas$marker == "Telomerase", 1:4]

# Simple visualization ...

plotdata(glas.t,                # Data frame
         two.by.two = FALSE     # Data is given as: (tp, n1, fp, n2)
         )

glas.m1 <- metadiag(glas.t,                # Data frame
                    two.by.two = FALSE,    # Data is given as: (tp, n1, fp, n2)
                    re = "normal",         # Random effects distribution
                    re.model = "DS",       # Random effects on D and S
                    link = "logit",        # Link function
                    sd.Fisher.rho   = 1.7, # Prior standard deviation of correlation
                    nr.burnin = 1000,      # Iterations for burnin
                    nr.iterations = 10000, # Total iterations
                    nr.chains = 2,         # Number of chains
                    r2jags = TRUE)         # Use r2jags as interface to jags


 summary(glas.m1, digit=3)

 plot(glas.m1,                    # Fitted model
      level = c(0.5, 0.75, 0.95), # Credibility levels
      parametric.smooth = TRUE)   # Parametric curve


# Plot results: based on a non-parametric smoother of the posterior predictive rates .......

plot(glas.m1,                    # Fitted model
     level = c(0.5, 0.75, 0.95), # Credibility levels
     parametric.smooth = FALSE)  # Non-parametric curve


# Using the pipe command in the package dplyr ...............................................

library(dplyr)

glas.t %>%
  metadiag(re = "normal", re.model ="SeSp") %>%
  plot(parametric.smooth = FALSE, color.pred.points = "red")



# Visualization of posteriors of hyper-parameters .........................................
library(ggplot2)
library(GGally)
library(R2jags)
attach.jags(glas.m1)
hyper.post <- data.frame(mu.D, mu.S, sigma.D, sigma.S, rho)
ggpairs(hyper.post,                  # Data frame
       title = "Hyper-Posteriors",          # title of the graph
       lower = list(continuous = "density") # contour plots
       )


#............................................................................

# List of different statistical models:
#    1) Different link functions: logit, cloglog and probit

#    2) Different parametrization of random effects in the link scale:
#         DS = "differences of TPR and FPR"
#         SeSp = "Sensitivity and Specificity"

#    3) Different random effects distributions:
#       "normal" or "sm = scale mixtures".

#    4) For the scale mixture random effects:
#       split.w = TRUE => "split the weights".

#    5) For the scale mixture random effects:
#       df.estimate = TRUE => "estimate the degrees of freedom".

#    6) For the scale mixture random effects:
#       df.estimate = TRUE => "estimate the degrees of freedom".

#    7) For the scale mixture random effects:
#       df = 4 => "fix the degrees of freedom to a particual value".
#       Note that df = 1 fits a Cauchy bivariate distribution to the random effects.

# logit-normal-DS
m <- metadiag(glas.t, re = "normal", re.model = "DS", link = "logit")
summary(m)
plot(m)

# cloglog-normal-DS
summary(metadiag(glas.t, re = "normal", re.model = "DS", link = "cloglog"))

# probit-normal-DS
summary(metadiag(glas.t, re = "normal", re.model = "DS", link = "probit"))
# logit-normal-SeSp
summary(metadiag(glas.t, re = "normal", re.model = "SeSp", link = "logit"))

# cloglog-normal-SeSp
summary(metadiag(glas.t, re = "normal", re.model = "SeSp", link = "cloglog"))
# probit-normal-SeSp
summary(metadiag(glas.t, re = "normal", re.model = "SeSp", link = "probit"))

# logit-sm-DS
summary(metadiag(glas.t, re = "sm", re.model = "DS", link = "logit", df = 1))

# cloglog-sm-DS
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "cloglog", df = 1))
plot(m, parametric.smooth = FALSE)

# probit-sm-DS
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "probit", df = 1))
plot(m, parametric.smooth = FALSE)

# logit-sm-SeSp
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "logit", df = 1))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# cloglog-sm-SeSp
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "cloglog", df = 1))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# probit-sm-SeSp
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "probit", df = 1))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# logit-sm-DS-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "logit",
 df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# cloglog-sm-DS-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "cloglog",
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# probit-sm-DS-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "probit",
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# logit-sm-SeSp-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "probit",
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# cloglog-sm-SeSp-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "cloglog",
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# probit-sm-SeSp-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "probit",
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))

# split.w ...................................................................

# logit-sm-DS
summary(m <- metadiag(glas.t, re = "sm", re.model = "DS", link = "logit", split.w = TRUE, df = 10))
plot(m)

# cloglog-sm-DS
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "cloglog", split.w = TRUE, df = 4))
plot(m)

# probit-sm-DS
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "probit", split.w = TRUE, df = 4))
plot(m, parametric.smooth = FALSE)

# logit-sm-SeSp
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "logit", split.w = TRUE, df = 1))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

# cloglog-sm-SeSp
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "cloglog", split.w = TRUE, df = 1))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

# probit-sm-SeSp
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "probit", split.w = TRUE, df = 1))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)


# logit-sm-DS-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "logit", split.w = TRUE,
 df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

# cloglog-sm-DS-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "cloglog", split.w = TRUE,
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

# probit-sm-DS-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "DS", link = "probit", split.w = TRUE,
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

# logit-sm-SeSp-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "probit", split.w = TRUE,
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

# cloglog-sm-SeSp-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "cloglog", split.w = TRUE,
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

# probit-sm-SeSp-df
summary(m<-metadiag(glas.t, re = "sm", re.model = "SeSp", link = "probit", split.w = TRUE,
df.estimate = TRUE))
plot(m, parametric.smooth = FALSE, level = c(0.5, 0.9))
plotw(m)

}