# Bayesian Meta-Analysis of diagnostic test data 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**. ```r 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 ) ``` ## 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. ## Returns 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. ## 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. ## Examples ```r ## Not run: # 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) ## End(Not run) ``` ## 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.

metadiag function