```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

library(BMA)

library(metafor)

library(mada)

library(ggplot2)

1 Load dataset

ebpp <- read.csv(“C:/DATA/EBPP.csv”, header = TRUE, fileEncoding = “latin1”)

2 Apply continuity correction

ebpp\(TP_adj <- ebpp\)TP + 0.5

ebpp\(FN_adj <- ebpp\)FN + 0.5

ebpp\(TN_adj <- ebpp\)TN + 0.5

ebpp\(FP_adj <- ebpp\)FP + 0.5

3 Calculate Se and Sp

ebpp\(Se <- ebpp\)TP_adj / (ebpp\(TP_adj + ebpp\)FN_adj)

ebpp\(Sp <- ebpp\)TN_adj / (ebpp\(TN_adj + ebpp\)FP_adj)

4 Logit transform

ebpp\(logitSe <- qlogis(ebpp\)Se)

ebpp\(logitSp <- qlogis(ebpp\)Sp)

5 Variance of logit

ebpp\(varLogitSe <- 1 / ebpp\)TP_adj + 1 / ebpp$FN_adj

ebpp$varLogitSp <- 1 / ebpp\(TN_adj + 1 / ebpp\)FP_adj

res_se <- rma(yi = logitSe, vi = varLogitSe, data = ebpp, method = “REML”)

res_sp <- rma(yi = logitSp, vi = varLogitSp, data = ebpp, method = “REML”)

se_summary <- plogis(coef(res_se))

sp_summary <- plogis(coef(res_sp))

se_CI <- plogis(confint(res_se))

p_CI <- plogis(confint(res_sp))

data.frame( Parameter = c(“Sensitivity”, “Specificity”), Estimate = c(se_summary, sp_summary), CI_Lower = c(se_CI[1], sp_CI[1]), CI_Upper = c(se_CI[2], sp_CI[2]), I2 = c(res_se\(I2, res_sp\)I2) )

forest(res_se, xlab = “Logit Sensitivity”, slab = ebpp$Study)

forest(res_sp, xlab = “Logit Specificity”, slab = ebpp$Study)

forest(res_se, xlab = “Sensitivity”, transf = plogis, slab = ebpp$Study)

forest(res_sp, xlab = “Specificity”, transf = plogis, slab = ebpp$Study)

ebpp\(LR_pos <- ebpp\)Se / (1 - ebpp$Sp)

ebpp$LR_neg <- (1 - ebpp$Se) / ebpp$Sp

ebpp\(logLR_pos <- log(ebpp\)LR_pos)

ebpp\(logLR_neg <- log(ebpp\)LR_neg)

ebpp\(varLogLR_pos <- 1 / ebpp\)TP_adj + 1 / ebpp$FP_adj

ebpp$varLogLR_neg <- 1 / ebpp$FN_adj + 1 / ebpp$TN_adj

res_lr_pos <- rma(yi = logLR_pos, vi = varLogLR_pos, data = ebpp, method = “REML”)

res_lr_neg <- rma(yi = logLR_neg, vi = varLogLR_neg, data = ebpp, method = “REML”)

LRp_summary <- exp(coef(res_lr_pos))

LRn_summary <- exp(coef(res_lr_neg))

AUC_approx <- (LRp_summary - LRn_summary + 1) / (2 * (LRp_summary - LRn_summary))

data.frame( Metric = c(“LR+”, “LR-”, “Approx AUC”), Estimate = round(c(LRp_summary, LRn_summary, AUC_approx), 3) )

plot(1 - ebpp\(Sp, ebpp\)Se, xlab = “1 - Specificity”, ylab = “Sensitivity”, xlim = c(0, 1), ylim = c(0, 1), pch = 16, col = “gray30”, main = “Summary ROC Plot”)

abline(0, 1, col = “gray”, lty = 2)

points(1 - sp_summary, se_summary, pch = 19, col = “red”, cex = 1.5)

legend(“bottomright”, legend = c(“Studies”, “Summary Estimate”), col = c(“gray30”, “red”), pch = c(16, 19))

mada_obj <- mada(TP = ebpp$TP, FP = ebpp$FP, FN = ebpp$FN, TN = ebpp$TN, names = ebpp$Study)

reitsma_model <- reitsma(mada_obj)

plot(reitsma_model, sroclwd = 2, main = “SROC Curve (Bivariate Model)”)

auc_biv <- AUC(reitsma_model) cat(“AUC (from bivariate model):”, round(auc_biv, 3), “”)