In JinjiPang/bstrp:

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "80%"
)

library(pROC)
library(tidyverse)
library(ggplot2)
library(boot)

1.Explore the data using pROC package

1.1 read in data

dfroc<-read.csv("../data/ROC.csv")

dfnew<-dfroc%>%
  mutate(
    y=ifelse(Exposed=="NEG", "0","1")
  )

1.2 calculate roc and visualization for all three assays

## ROC curves for wvELISA, IFA and SVN
rocall<-roc(Exposed~wvELISA+IFA+VNResult,data=dfroc)
plot1<-ggroc(rocall)
plot1

1.3 calculate AUC seperately using pROC package

## ROC curves for wvELISA
roc1<-roc(Exposed~wvELISA,data=dfnew)
## ROC curves for IFA 
roc2<-roc(Exposed~IFA,data=dfnew)
## ROC curves forSVN
roc3<-roc(Exposed~VNResult,data=dfnew)

##calculate AUC and 95% CI
CI<-c(ci.auc(roc1),ci.auc(roc2),ci.auc(roc3))
CI<-matrix(CI, nrow = 3, byrow=T)

AUC<-data.frame(CI[,2],CI[,1],CI[,3])
rownames(AUC)<-c("wvELISA","IFA","SVN")

AUC<-AUC%>%knitr::kable(caption="95% CI of AUC ", col.names = c("AUC","Lowwer","Upper"))
print(AUC)

1.4 pair-wise ROC test (test difference in AUC) using pROC package

## wv_ELISA vs IFA
roc.test(roc1, roc2, method="bootstrap", boot.n=10000)
## wv_ELISA vs SVN
roc.test(roc1, roc3, method="bootstrap", boot.n=10000)
## IFA vs SVN
roc.test(roc2, roc3, method="bootstrap", boot.n=10000)

2.Explore the data while considering the repeated measures data structure

2.1 calculate AUC seperatly using bootstrap methods (boot and pROC packages)

set.seed(0702)
## get the unique Pig ID
pig_id <- unique(dfnew$Pig.ID)

## function for compute the bootstrap statistics AUC
compute_auc <- function(data, indices, formula, data_all) {
  # id from bootstrap sampling
  selected_id <- data[indices]
  # 
  boot_data_list <- lapply(selected_id, function(id) {
    data_all[data_all$Pig.ID == id , ]
  })

  # get all data with selected id
  boot_data <- do.call(rbind, boot_data_list)

  return(roc(formula, data=boot_data)$auc)
}


#### compute AUC for wvELISA
result_elisa <- boot(data=pig_id, statistic=compute_auc,
                   R=2000, formula=y~wvELISA, data_all=dfnew)
#### compute AUC for  IFA
result_ifa <- boot(data=pig_id, statistic=compute_auc,
                       R=2000, formula=y~IFA, data_all=dfnew)
#### compute AUC for SVN
result_svn <- boot(data=pig_id, statistic=compute_auc,
                       R=2000, formula=y~VNResult, data_all=dfnew)



wvELISA<-quantile(result_elisa$t[,1], c(0.5,0.025, 0.975))
IFA<-quantile(result_ifa$t[,1], c(0.5, 0.025, 0.975))
SVN<-quantile(result_svn$t[,1], c(0.5, 0.025, 0.975))

df<-data.frame(wvELISA,IFA,SVN)
rownames(df)<-c("AUC","Lower","Upper")
print(df)

2.2. pair-wise test (test difference in AUC) using boot and pROC packages

## function for compute the bootstrap statistics AUC difference
compare_auc <- function(data, indices, data_all, formula1, formula2) {
  # id from bootstrap sampling
  selected_id <- data[indices]

  boot_data_list <- lapply(selected_id, function(id) {
    data_all[data_all$Pig.ID == id , ]
  })

  # get all data with selected id
  boot_data <- do.call(rbind, boot_data_list)

  auc1 <- roc(formula1, data=boot_data)$auc
  auc2 <- roc(formula2, data=boot_data)$auc

  return(auc1 - auc2)
}

#### compute AUC difference for wvELISA and IFA
auc_compare_elisa_ifa <- boot(data=pig_id, statistic=compare_auc,
                              R=2000, data_all= dfnew,
                              formula1=y~wvELISA, formula2 = y~IFA)

quantile(auc_compare_elisa_ifa$t[,1], c(0.025, 0.975))

diff1 <- auc_compare_elisa_ifa$t[,1]

ELISA_IFA_pvalue= (1 - pt( (mean(diff1) - 0)/sd(diff1), df = 1999)) * 2


#### compute AUC difference for wvELISA and SVN
auc_compare_elisa_svn <- boot(data=pig_id, statistic=compare_auc,
                              R=2000, data_all= dfnew,
                              formula1=y~wvELISA, formula2 = y~VNResult)

quantile(auc_compare_elisa_svn$t[,1], c(0.025, 0.975))

diff2 <- auc_compare_elisa_svn$t[,1]

ELISA_SVN_pvalue = (1 - pt( (mean(diff2) - 0)/ sd(diff2), df = 1999)) * 2


#### compute AUC difference for IFA and SVN
auc_compare_ifa_svn <- boot(data=pig_id, statistic=compare_auc,
                              R=2000, data_all= dfnew,
                              formula1=y~IFA, formula2 = y~VNResult)

quantile(auc_compare_ifa_svn$t[,1], c(0.025, 0.975))

diff3 <- auc_compare_ifa_svn$t[,1]

##note that on average, IFA's AUC < SVN's AUC
IFA_SVN_pvalue = (pt((mean(diff3) - 0)/ sd(diff3), df = 1999)) * 2

df<-data.frame(ELISA_IFA_pvalue, ELISA_SVN_pvalue,IFA_SVN_pvalue)

print(df)