rxTest({ if (!.Call(`_rxode2_isIntel`)) { test_that("test binomProb()", { x <- rbinom(1001, 1, prob = 0.05) m <- mean(x) v <- m * (1 - m) s <- sqrt(v) n <- 1001 z <- qnorm(1 - 0.025) z2 <- z * z c1 <- 2*n*m + z2 c2 <- z*sqrt(z2 - 1.0/n + 4*n*v + (4.0*m-2.0))+1.0 c3 <- 2*(n+z2) z2.5 <- (c1 + c(-1, 1) * c2) / c3 n <- 1001 z <- qnorm(1 - 0.05) z2 <- z * z c1 <- 2*n*m + z2 c2 <- z*sqrt(z2 - 1.0/n + 4*n*v + (4.0*m-2.0))+1.0 c3 <- 2*(n+z2) z5 <- (c1 + c(-1, 1) * c2) / c3 t1 <- c("2.5%"=z2.5[1], "5%"=z5[1], "50%"=m, "95%"=z5[2], "97.5%"=z2.5[2]) expect_equal(binomProbs(x, ciMethod="wilsonCorrect"), t1) expect_equal(binomProbs(x, ciMethod="wc"), t1) t2 <- c(c("mean"=m, "var"=v, "sd"=s, "n"=n), t1) expect_equal(binomProbs(x, onlyProbs=FALSE, ciMethod="wilsonCorrect"), t2) x2 <- c(x, NA_real_) setNames(rep(NA_real_, length(t1)),names(t1)) expect_equal(binomProbs(x2, ciMethod="wilsonCorrect"), setNames(rep(NA_real_, length(t1)),names(t1))) expect_equal(binomProbs(x2, onlyProbs=FALSE, ciMethod="wilsonCorrect"), setNames(rep(NA_real_, length(t2)),names(t2))) expect_equal(binomProbs(x2,na.rm=TRUE, ciMethod="wilsonCorrect"), t1) expect_equal(binomProbs(x2, onlyProbs=FALSE, na.rm=TRUE, ciMethod="wilsonCorrect"), t2) n <- 50 z <- qnorm(1 - 0.025) z2 <- z * z c1 <- 2*n*m + z2 c2 <- z*sqrt(z2 - 1.0/n + 4*n*v + (4.0*m-2.0))+1.0 c3 <- 2*(n+z2) z2.5 <- (c1 + c(-1, 1) * c2) / c3 n <- 50 z <- qnorm(1 - 0.05) z2 <- z * z c1 <- 2*n*m + z2 c2 <- z*sqrt(z2 - 1.0/n + 4*n*v + (4.0*m-2.0))+1.0 c3 <- 2*(n+z2) z5 <- (c1 + c(-1, 1) * c2) / c3 t1 <- c("2.5%"=z2.5[1], "5%"=z5[1], "50%"=m, "95%"=z5[2], "97.5%"=z2.5[2]) expect_equal(binomProbs(x, n=50, ciMethod="wilsonCorrect"), t1) # now wilson x <- rbinom(1001, 1, prob = 0.05) getVals <- function(x, a=0.025) { p <- mean(x) n <- length(x) z <- qnorm(1 - a) z2 <- z ^ 2 c0 <- 1.0 / (1.0 + z2 / n) c1 <- c0 * (p + z2 / (2 * n)) c2 <- z * c0 * sqrt(p * (1 - p) / n + z2 / (4 * n * n)) (c1 + c(-1, 1) * c2) } m <- mean(x) z2.5 <- getVals(x, 0.025) z5 <- getVals(x, 0.05) t1 <- c("2.5%"=z2.5[1], "5%"=z5[1], "50%"=m, "95%"=z5[2], "97.5%"=z2.5[2]) expect_equal(binomProbs(x, ciMethod="wilson"), t1) # wald x <- rbinom(1001, 1, prob = 0.05) getVals <- function(x, a=0.025) { p <- mean(x) n <- length(x) z <- qnorm(1 - a) c1 <- p c2 <- z * sqrt(p * (1 - p) / n) (c1 + c(-1, 1) * c2) } m <- mean(x) z2.5 <- getVals(x, 0.025) z5 <- getVals(x, 0.05) t1 <- c("2.5%"=z2.5[1], "5%"=z5[1], "50%"=m, "95%"=z5[2], "97.5%"=z2.5[2]) expect_equal(binomProbs(x, ciMethod="wald"), t1) # Agresti-Coull x <- rbinom(1001, 1, prob = 0.05) getVals <- function(x, a=0.025) { p <- mean(x) n <- length(x) z <- qnorm(1 - a) z2 <- z ^ 2 nh <- n + z2 ns <- p * n ph <- 1 / nh * (ns + z2 / 2) c1 <- ph c2 <- z * sqrt(ph * (1 - ph) / nh) (c1 + c(-1, 1) * c2) } m <- mean(x) z2.5 <- getVals(x, 0.025) z5 <- getVals(x, 0.05) t1 <- c("2.5%"=z2.5[1], "5%"=z5[1], "50%"=m, "95%"=z5[2], "97.5%"=z2.5[2]) expect_equal(binomProbs(x, ciMethod="agrestiCoull"), t1) expect_equal(binomProbs(x, ciMethod="ac"), t1) expect_equal(names(binomProbs(x, pred=TRUE)), names(binomProbs(x))) }) } })