# Replication script for the examples of Section 4 from: # "bamdit: An R Package for Bayesian Meta-Analysis of Diagnostic Test Data" # Begin of the script ... set.seed(2017) ## ------------------------------------------------------------------------ library("bamdit") data("glas") head(glas) plotdata(glas, # Data frame group = glas$marker, # groupping variable max.size = 20) # scale of circles ## ------------------------------------------------------------------------ glas.t <- glas[glas$marker == "Telomerase", 1:4] ## ------------------------------------------------------------------------ plotdata(glas.t) ## ------------------------------------------------------------------------ glas.m1 <- metadiag(glas.t, # Data frame 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 = 4, # Number of chains r2jags = TRUE) # Use r2jags as interface to jags ## ------------------------------------------------------------------------ summary(glas.m1, digits = 3) ## ------------------------------------------------------------------------ library("R2jags") attach.jags(glas.m1) cor(se.pool, sp.pool) ## ------------------------------------------------------------------------ plot(glas.m1, # Fitted model level = c(0.5, 0.75, 0.95), # Credibility levels parametric.smooth = TRUE) # Parametric curve ## ------------------------------------------------------------------------ plotsesp(glas.m1) ## ------------------------------------------------------------------------ bsroc(glas.m1, # Fitted model level = c(0.025, 0.5, 0.975), # Credibility levels plot.post.bauc = TRUE, # Include the posterior of the AUC fpr.x = seq(0.01, 0.75, 0.01), # Grid of values for FPR lower.auc = 0.01, # Lower limit for the BAUC upper.auc = 0.75, # Upper limit for the BAUC partial.AUC = FALSE) ## ------------------------------------------------------------------------ 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 with MCMC realizations title = "Hyper-Posteriors", # title of the graph lower = list(continuous = "density") # contour plots ) ## ------------------------------------------------------------------------ glas.m2 <- metadiag(glas.t, # Data frame re = "sm", # Scale mixture of normals link = "logit", # Link function sd.Fisher.rho = 1.7, # Prior standard deviation of correlation df.estimate = TRUE, # Degrees of freedom estimated from the data split.w = TRUE, # Different weights for each component nr.burnin = 10000, # Iterations for burnin nr.iterations = 100000, # Total iterations nr.chains = 1, # Number of chains r2jags = TRUE) # Use r2jags as interface to jags ## ------------------------------------------------------------------------ glas.m2 ## ------------------------------------------------------------------------ plotw(m = glas.m2) ## ------------------------------------------------------------------------ glas.t[c(5, 7, 10), ] ## ------------------------------------------------------------------------ dat.hat <- data.frame(tpr = glas.t[ , 1]/glas.t[ , 2], fpr = glas.t[ , 3]/glas.t[, 4], n = glas.t[ , 2] + glas.t[ , 4]) dat.hat[c(5, 7, 10), ] ## ------------------------------------------------------------------------ plotcompare(m1 = glas.m1, # Model 1 object m2 = glas.m2, # Model 2 object m1.name = "Binomial + Normal", # Label for Model 1 m2.name = "Binomial + Scale mixtures", # Label for Model 2 level = 0.95) ## ------------------------------------------------------------------------ data("ct") gr <- with(ct, factor(design, labels = c("Retrospective study", "Prospective study"))) plotdata(ct, # Data frame group = gr, # Groupping variable y.lo = 0.75, # Lower limit of y-axis x.up = 0.75, # Upper limit of x-axis alpha.p = 0.5, # Transparency of the balls max.size = 20) # Scale the circles ## ------------------------------------------------------------------------ ct.m <- metadiag(ct, re = "sm", # Scale mixture of normals link = "logit", # Link function df.estimate = TRUE, # Degrees of freedom estimated from the data split.w = TRUE, # Different weights for each component nr.burnin = 1000, # Iterations for burnin nr.iterations = 10000, # Total iterations nr.chains = 4, # Number of chains r2jags = TRUE) # Use r2jags as interface to jags ## ------------------------------------------------------------------------ plotw(m = ct.m, # The fitted model group = gr # The groupping factor ) ## ------------------------------------------------------------------------ m1.ct <- metadiag(ct[ct$design==1, 1:4]) # Restrospective studies m2.ct <- metadiag(ct[ct$design==2, 1:4]) # Prospective studies plotcompare(m1.ct, m2.ct, m1.name = "Retrospective design", m2.name = "Prospective design", group = gr, limits.x = c(0, 0.75), limits.y = c(0.65, 1)) # End of the script ...