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]))) }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.