ns <- matrix(NA, 5, 2)
ns[,1] <- c(20, 40, 60, 80, 100)
ns[,2] <- c(100, 80, 60, 40, 20)
ps <- c(0.05, 0.2, 0.5, 0.8, 0.95)
ovresult <- NULL

set.seed(1234)
nrep <- 1000
for(k in 1:nrow(ns)) {
	n <- ns[k,]
	for(j in 1:length(ps)) {
		p <- ps[j]
		result <- rep(NA, nrep)
		for(i in 1:nrep) {
			y1 <- rbinom(n[1], 1, p)
			y2 <- rbinom(n[2], 1, p)
			yes1 <- sum(y1)
			yes2 <- sum(y2)
			total1 <- length(y1)
			total2 <- length(y2)
			out <- prop.test(x = c(yes1, yes2), n = c(total1, total2))
			result[i] <- out$p.value
		}
		result <- cbind(k, j, result)
		ovresult <- rbind(ovresult, result)
	}
}

colnames(ovresult) <- c("n", "prop", "p")
ovresult <- data.frame(ovresult, sig = ovresult[,"p"] < 0.05)
ovresult[,"n"] <- factor(ovresult[,"n"], labels = c("20,100", "40,80", "60,60", "80,40", "100,20"))
ovresult[,"prop"] <- factor(ovresult[,"prop"], labels = c("0.05", "0.2", "0.5", "0.8", "0.95"))
aggregate(sig ~ n + prop, data = ovresult, FUN = mean)