In vqf/mtest: Direct statistical test to compare multinomial trials

library(mtest)
library(Barnard)
set.seed(42)
MAKEPDF <- T
PDFPATH <- 'D:/Math/Bernouilli/ms/vtex-soft-texsupport.ims-aos-07708a9/vtex-soft-texsupport.ims-aos-07708a9/fig2.pdf'
texp <- list(c(10, 0), c(7, 0))
nmc <- 200000
mt <- mcm(texp, nmc, algo=6, NRUNIF = 0) #Alternative model, positives are true
mf <- mcm(texp, nmc, algo=4, NRUNIF = 0) #Null model, positives are false
mf <- calCumul(mf) 
mt <- calCumul(mt)
tf <- bernDist(texp)
tf <- calCumul(tf)
names(tf) <- mapply(function(x){tmp <- paste(as.vector(x$desc), collapse = '_'); names(x) <- tmp}, tf)
## Add Barnard p-vals
for (nm in names(tf)){
  mtx <- tf[[nm]]$desc
  btest <- Barnard::barnard.test(mtx[1, 1], mtx[2, 1], mtx[1, 2], mtx[2, 2])
  ftest <- as.numeric(fisher.test(mtx)$p.value)
  pval <- as.numeric(btest$p.value[[2]])
  tf[[nm]]$bpval <- pval
  tf[[nm]]$fpval <- ftest
}
tr <- matrix(ncol = 2, nrow = 0)
for (n in names(mf)){
  tr <- rbind(tr, c(tf[[n]]$cval, mf[[n]]$cval ))
}
colnames(tr) <- c('Theoretical', 'Simulated')
smoothScatter(tr)
title('2x2 m-test p-value, null hypothesis')
abline(0, 1, col='red', lty=2)

tr <- matrix(ncol = 2, nrow = 0)
for (n in names(mt)){
  tr <- rbind(tr, c(tf[[n]]$cval, mt[[n]]$cval ))
}
colnames(tr) <- c('Theoretical', 'Simulated')
smoothScatter(tr)
title('2x2 m-test p-value, alternative hypothesis')
abline(0, 1, col='red', lty=2)

rtab <- matrix(ncol = 7, nrow = 0)
flt <- function(th, sm, alpha, tp, entry="cval"){
  result <- 0
  total <- 0
  all <- names(sm)
  for (x in all){
    success <- F
    val <- as.numeric(th[[x]][entry])
    if ((tp == 'fp' || tp == 'tp') && val < a){
      success <- T
    }
    if ((tp == 'fn' || tp == 'tn') && val >= a){
      success <- T
    }
    if (success == T) result <- result + sm[[x]]$n
    total <- total + sm[[x]]$n
  }
  result <- result / total
  return(result)
}

colnames(rtab) <- c('alpha', 'fp_mtest', 'tp_mtest', 
                    'fp_barnard', 'tp_barnard',
                    'fp_fisher', 'tp_fisher')
for (a in seq(0, 1, 0.005)){
  mfp <- flt(tf, mf, a, 'fp')
  mtp <- flt(tf, mt, a, 'tp')
  bfp <- flt(tf, mf, a, 'fp', 'bpval')
  btp <- flt(tf, mt, a, 'tp', 'bpval')
  ffp <- flt(tf, mf, a, 'fp', 'fpval')
  ftp <- flt(tf, mt, a, 'tp', 'fpval')

  rtab <- rbind(rtab, c(a, mfp, mtp, bfp, btp, ffp, ftp))

}
par(pin=c(3, 3))
plot(rtab[,2], rtab[,3], xlim=c(0,1), ylim=c(0,1), type='l', xlab = 'FPR', ylab='TPR', col="black")
title('ROC')
lines(rtab[,4], rtab[,5], lty=2, col="red")
lines(rtab[,6], rtab[,7], lty=3, col="blue")
legend(legend = c("m-test", "Barnard", "Fisher"), lty = c(1, 2, 2), col=c('black', 'red', 'blue'), x='right', y='middle')
plot(rtab[,1], rtab[,3], xlim=c(0,1), ylim=c(0,1), type='l', xlab = 'alpha', ylab='TPR', col="black")
title('Power')
lines(rtab[,1], rtab[,5], lty=2, col="red")
lines(rtab[,1], rtab[,7], lty=3, col="blue")
legend(legend = c("m-test", "Barnard", "Fisher"), lty = c(1, 2, 2), col=c('black', 'red', 'blue'), x='right', y='middle')

if (MAKEPDF == T){
  pdf(PDFPATH, width = 9, height = 5)
  par(mfrow=c(1, 2), pty="s", oma=c(0, 1.1, 0, 1), pin=c(2.5, 2.5))
  plot(rtab[,2], rtab[,3], xlim=c(0,1), ylim=c(0,1), type='l', xlab = 'FPR', ylab='TPR', col="black")
  title('ROC')
  lines(rtab[,4], rtab[,5], lty=2, col="red")
  lines(rtab[,6], rtab[,7], lty=2, col="blue")
  legend(legend = c("m-test", "Barnard", "Fisher"), lty = c(1, 2, 2), col=c('black', 'red', 'blue'), x='right', y='middle')
  par(pin=c(2.5, 2.5))
  plot(rtab[,1], rtab[,3], xlim=c(0,1), ylim=c(0,1), type='l', xlab = 'alpha', ylab='TPR', col="black")
  title('Power')
  lines(rtab[,1], rtab[,5], lty=2, col="red")
  lines(rtab[,1], rtab[,7], lty=2, col="blue")
  legend(legend = c("m-test", "Barnard", "Fisher"), lty = c(1, 2, 2), col=c('black', 'red', 'blue'), x='right', y='middle')
  dev.off()
}