knitr::opts_chunk$se
## function (...) 
## {
##     dots = resolve(...)
##     if (length(dots)) 
##         defaults <<- merge(dots)
##     invisible(NULL)
## }
## <bytecode: 0x00000000138edf70>
## <environment: 0x0000000013903210>
## Code for BBN meta-analysis and Prediction region
library(readxl)
library(lme4)
## Warning: package 'lme4' was built under R version 4.0.5
## Loading required package: Matrix
library(mada)
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 4.0.5
## Loading required package: ellipse
## 
## Attaching package: 'ellipse'
## The following object is masked from 'package:graphics':
## 
##     pairs
## Loading required package: mvmeta
## This is mvmeta 1.0.3. For an overview type: help('mvmeta-package').
Pancreas <- read_excel("liquidPancreas.xlsx")

X <- Pancreas
X<-as.data.frame(X)
N <- length(X$TP)

# Set up the data
# Generate 5 new variables of type long. We need these before we can reshape the data.
# n1 is number diseased
# n0 is number without disease
# true1 is number of true positives
# true0 is the number of true negatives
# study is the unique identifier for each study. 

X$n1 <- X$TP+X$FN
X$n0 <- X$FP+X$TN
X$true1 <- X$TP
X$true0 <- X$TN 

X$study <- 1:N

# Reshape the data from wide to long format 

Y = reshape(X, direction = "long", varying = list( c("n1" , "n0") , c( "true1","true0" ) ) ,
            timevar = "sens" , times = c(1,0) , v.names = c("n","true") ) 

# Sort data by study to cluster the 2 records per study together 
Y = Y[order(Y$id),]  
Y$spec<- 1-Y$sens

Perform meta-analysis # Perform meta-analysis

MA_Y = glmer( formula = cbind(  true , n - true ) ~ 0 + sens + spec + (0+sens + spec|study),
              data = Y , family = binomial , nAGQ=1 , verbose=0 )
             
 # More detail can be obtained by using the summary command 
ma_Y = summary(MA_Y)
labels(ma_Y) 
##  [1] "methTitle"    "objClass"     "devcomp"      "isLmer"       "useScale"    
##  [6] "logLik"       "family"       "link"         "ngrps"        "coefficients"
## [11] "sigma"        "vcov"         "varcor"       "AICtab"       "call"        
## [16] "residuals"    "fitMsgs"      "optinfo"
ma_Y$coeff 
##      Estimate Std. Error  z value     Pr(>|z|)
## sens 1.552911  0.3376913 4.598611 4.253178e-06
## spec 2.365975  0.2766763 8.551421 1.215826e-17
# Logit sensitivity and specificity
lsens = ma_Y$coeff[1,1]
lspec = ma_Y$coeff[2,1]
# Calculate confidence and predictive regions (based on Harbord et al 2005)

# Start by calculating derived and set parameters
seB     <- ma_Y$coefficients[2,2]
seA     <- ma_Y$coefficients[1,2]
r       <- ma_Y$vcov@x[2] / (seA*seB)
varA    <- ma_Y$varcor$study[1,1]
varB    <- ma_Y$varcor$study[2,2]
sAB     <- ma_Y$varcor$study[1,2]
covAB   <- ma_Y$vcov@x[2]
sepredA <- sqrt(varA + seA^2)
sepredB <- sqrt(varB + seB^2)
rpredAB <- (sAB + covAB) / (sepredA*sepredB)
level   <- 95
f       <- qf(0.95, df1=2, df2=N-2)
croot   <- sqrt(2*f)


# Empty data frames
conf_region <- (rep(NA, 361))
pred_region <- (rep(NA, 361))

# Confidence region
for (i in seq(0, 2*pi, length.out=361)){
  confB <- lspec + (seB*croot*cos(i))
  confA <- lsens + (seA*croot*cos(i + acos(r)))
  confsens <- plogis(confA)
  confspec <- plogis(confB)
  conf_i <- data.frame(X=1-confspec, Y=confsens)
  # Add results for most recent value to a data frame which contains the results of previous values
  conf_region<-rbind(conf_region, conf_i)
}
conf_region <- conf_region[2:362,]

# Predictive region
for (i in seq(0, 2*pi, length.out=361)){
  predB <- lspec + (sepredB*croot*cos(i))
  predA <- lsens + (sepredA*croot*cos(i + acos(rpredAB)))
  predsens <- plogis(predA)
  predspec <- plogis(predB)
  pred_i <- data.frame(X=1-predspec, Y=predsens)
  # Add results for most recent value to a data frame which contains the results of previous values
  pred_region<-rbind(pred_region, pred_i)
}
pred_region <- pred_region[2:362,]

# Sensitivity and specificity calculations for each study
X$sens <- X$TP / (X$TP + X$FN)
X$spec <- X$TN / (X$FP + X$TN)


study_level <- data.frame(ID=X$study, TP=X$TP, FN=X$FN, FP=X$FP, TN=X$TN, N=(X$TP+X$FN+X$FP+X$TN)
                          ,Sensitivity=X$sens, Specificity=X$spec)

# Mean point
Sens = plogis(lsens) 
Spec = plogis(lspec) 


mean_point=data.frame(x=Spec, y=Sens)

Including Plots

SROC plot with prediction and confidence regions

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
c<-ggplot(study_level,aes(Specificity,Sensitivity)) + geom_path(aes((1-X),Y,lty="95% prediction region"),data=pred_region)+
     geom_path(aes((1-X),Y,lty="95% confidence region"), data=conf_region)
c + geom_point() +xlim(1,0)+ylim(0,1)+geom_point(aes(x,y,col="Summary Estimate"),data=mean_point)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.