library(mtest)
calgo <- 4
set.seed(42)
ex <- list(c(17, 0), c(11, 0))
# MC Simulation
m <- mcm(ex, 100000, algo=calgo, NRUNIF = 0)
m <- calCumul(m)
# Theoretical m-test
f <- bernDist(ex)
f <- calCumul(f)
names(f) <- mapply(function(x){tmp <- paste(as.vector(x$desc), collapse = '_'); names(x) <- tmp}, f)
# Side by side p-values
r <- matrix(nrow = 0, ncol = 2)
for (n in names(m)){
  if (!is.null(f[[n]]) && !is.null(m[[n]])){
    r <- rbind(r, c(f[[n]]$cval, m[[n]]$cval))
  }
}
k <- lm(r[,2] ~ r[,1] + 0)
k

colnames(r) <- c('Theoretical', 'Simulated')
smoothScatter(r)
title('m-test p-value')
abline(0, 1, col='red', lty=2)
abline(k, col="black", lty=1)
legend(legend = c("Regression", "Slope=1"), lty = c(1, 2), col=c('black', 'red'), x='right', y='middle')
library(Barnard)
n1 <- 17
n2 <- 11
exb <- list(c(n1, 0), c(n2, 0))
# MC Simulation
mb <- mcm(exb, 100000, algo=calgo, NRUNIF = 0)
mb <- calCumul(mb)
s <- matrix(nrow = 0, ncol = 2)
for (i in 1:n1){
  for (j in 1:n2){
    dsc <- c(i, j, n1-i, n2-j)
    nm <- paste(dsc, collapse = '_')
    br <- barnard.test(dsc[1], dsc[2], dsc[3], dsc[4])
    if (!is.null(br$p.value[2]) && !is.null(m[[n]])){
      v <- as.numeric(br$p.value[[2]])
      w <- mb[[nm]]$cval
      s <- rbind(s, c(v, w))
    }
  }
}
colnames(s) <- c('Theoretical', 'Simulated')
smoothScatter(s)
title('Barnard test p-value')
kb <- lm(s[,2] ~ s[,1] + 0)
kb
abline(0, 1, col='red', lty=2)
abline(kb, col="black", lty=1)
legend(legend = c("Regression", "Slope=1"), lty = c(1, 2), col=c('black', 'red'), x='right', y='middle')
mnc <- ncol(m[[1]]$desc)
toh <- Reduce(init = matrix(ncol = mnc, nrow = 0), function(x, y){rbind(x, matrix(rep(as.vector(y[[1]]), y[[2]]), ncol = mnc, byrow = F))}, Map(function(x){list(as.matrix(x$desc), x$n)}, m))
for (i in 1:ncol(toh)){
  hist(toh[,i], breaks = c(-1:max(toh[,i])))
}


vqf/mtest documentation built on Dec. 23, 2021, 4:11 p.m.