tests/testthat/test-binomProb.R

In nlmixr2/rxode2: Facilities for Simulating from ODE-Based Models

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)))

    })
  }
})