MVP?

Author: pearsonca

DESCRIPTION

@@ -11,3 +11,6 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
+Suggests: 
+    testthat (>= 3.0.0)
+Config/testthat/edition: 3

NAMESPACE

@@ -1,2 +1,4 @@
 # Generated by roxygen2: do not edit by hand
 
+export(vfisher.test)
+import(data.table)

R/vfisher.R

@@ -26,31 +26,16 @@ pnhyper <- function(
   q, ncp = 1, upper.tail = FALSE, m, n, k, lo, hi, support, logdc
 ) {
 
-  lims <- which(ncp == 0 | is.infinite(ncp))
-  ones <- which(ncp == 1)
-  ncp1 <- which(ncp == 1)
-  res <- numeric(length(q))
-
-  res[lims] <- ifelse(
-    ncp == 0,
-    as.numeric(if (upper.tail) q <= lo else q >= lo),
-    as.numeric(if (upper.tail) q <= hi else q >= hi)
-  )
-
-  res[ones] <- phyper(
-    x[ones] - upper.tail, m[ones], n[ones], k[ones],
-    lower.tail = !upper.tail
-  )
-
-  res[-c(ones, lims)] <- mapply(
-    function(q, ncp, logdc, support) {
-      sum(dnhyper(ncp, logdc, support)[if (upper.tail) support >= q else support <= q])
-    },
-    q = q[-c(ones, lims)], ncp = ncp[-c(ones, lims)], logdc = logdc[-c(ones, lims)],
-    support = support[-c(ones, lims)]
-  )
+  if (ncp == 0) {
+    return(as.numeric(if (upper.tail) q <= lo else q >= lo))
+  } else if (is.infinite(ncp)) {
+    return(as.numeric(if (upper.tail) q <= hi else q >= hi))
+  } else if (ncp == 1) {
+    return(phyper(q - upper.tail, m, n, k, lower.tail = !upper.tail))
+  } else {
+    return(sum(dnhyper(ncp, logdc, support)[if (upper.tail) q <= support else support <= q]))
+  }
 
-  return(res)
 }
 
 # ncp.U == ncp_ci(..., lower = FALSE)
@@ -126,6 +111,8 @@ vfisher.test <- function(
     or <- rep(or, length.out = length(a))
   }
 
+  alternative <- match.arg(alternative)
+
   # matrix =
   # a, b
   # c, d
@@ -151,6 +138,20 @@ vfisher.test <- function(
     logdc :=  .(list(dhyper(support[[1]], m, n, k, log = TRUE))), by = rowid
   ]
 
+  if (alternative %in% c("less", "greater")) {
+    result_dt[, p.value := pnhyper(
+      a, or, upper.tail = alternative == "greater", m, n, k, lo, hi, support[[1]], logdc[[1]]
+    ), by = rowid]
+  } else {
+    result_dt[or == 0, p.value := as.numeric(a == lo)]
+    result_dt[is.infinite(or), p.value := as.numeric(a == hi)]
+    result_dt[!(or == 0 | is.infinite(or)), p.value := {
+      relErr <- 1 + 10^(-7)
+      d <- dnhyper(or, logdc[[1]], support[[1]])
+      sum(d[d <= d[a - lo + 1] * relErr])
+    }, by = rowid]
+  }
+
   # x == a
 
   result_dt[a == lo, estimate := 0]
@@ -164,16 +165,12 @@ vfisher.test <- function(
 
   if (conf.int) {
 
-    alternative <- match.arg(alternative)
-
     if (!((length(conf.level) == 1L) && is.finite(conf.level) &&
           (conf.level > 0) && (conf.level < 1)))
       stop("'conf.level' must be a single number between 0 and 1")
 
     sdcols <- c(sdcols, c("ci.lo", "ci.hi"))
 
-    setattr(result_dt, "conf.level", conf.level)
-
     if (alternative == "less") {
       result_dt[, ci.lo := 0]
       result_dt[, ci.hi := ncp_ci(a, 1 - conf.level, m, n, k, lo, hi, support[[1]], logdc[[1]], lower = FALSE), by = rowid]
@@ -188,68 +185,10 @@ vfisher.test <- function(
 
   }
 
-  result_dt[, .SD, .SDcols = sdcols]
-
-}
+  res <- result_dt[, .SD, .SDcols = sdcols]
 
+  if (conf.int) setattr(result_dt, "conf.level", conf.level)
 
-function (x, y = NULL, workspace = 2e+05, hybrid = FALSE, hybridPars = c(expect = 5,
-                                                                         percent = 80, Emin = 1), control = list(), or = 1, alternative = "two.sided",
-          conf.int = TRUE, conf.level = 0.95, simulate.p.value = FALSE,
-          B = 2000)
-{
-
-  PVAL <- NULL
-  else {
-    PVAL <- switch(alternative, less = pnhyper(x, or),
-                   greater = pnhyper(x, or, upper.tail = TRUE),
-                   two.sided = {
-                     if (or == 0) as.numeric(x == lo) else if (or ==
-                                                               Inf) as.numeric(x == hi) else {
-                                                                 relErr <- 1 + 10^(-7)
-                                                                 d <- dnhyper(or)
-                                                                 sum(d[d <= d[x - lo + 1] * relErr])
-                                                               }
-                   })
-  }
+  return(res)
 
-  ESTIMATE <- c(`odds ratio` = mle(x))
-  if (conf.int) {
-    ncp.U <- function(x, alpha) {
-      if (x == hi)
-        return(Inf)
-      p <- pnhyper(x, 1)
-      if (p < alpha)
-        uniroot(function(t) pnhyper(x, t) - alpha,
-                c(0, 1))$root
-      else if (p > alpha)
-        1/uniroot(function(t) pnhyper(x, 1/t) - alpha,
-                  c(.Machine$double.eps, 1))$root
-      else 1
-    }
-    ncp.L <- function(x, alpha) {
-      if (x == lo)
-        return(0)
-      p <- pnhyper(x, 1, upper.tail = TRUE)
-      if (p > alpha)
-        uniroot(function(t) pnhyper(x, t, upper.tail = TRUE) -
-                  alpha, c(0, 1))$root
-      else if (p < alpha)
-        1/uniroot(function(t) pnhyper(x, 1/t, upper.tail = TRUE) -
-                    alpha, c(.Machine$double.eps, 1))$root
-      else 1
-    }
-    CINT <- switch(alternative, less = c(0, ncp.U(x,
-                                                  1 - conf.level)), greater = c(ncp.L(x, 1 - conf.level),
-                                                                                Inf), two.sided = {
-                                                                                  alpha <- (1 - conf.level)/2
-                                                                                  c(ncp.L(x, alpha), ncp.U(x, alpha))
-                                                                                })
-    attr(CINT, "conf.level") <- conf.level
-  }
-  RVAL <- c(RVAL, list(conf.int = if (conf.int) CINT, estimate = ESTIMATE,
-                       null.value = NVAL))
-}
-structure(c(RVAL, alternative = alternative, method = METHOD,
-            data.name = DNAME), class = "htest")
 }

man/vfisher.test.Rd

@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview
+# * https://testthat.r-lib.org/articles/special-files.html
+
+library(testthat)
+library(vfisher)
+
+test_check("vfisher")

tests/testthat.R

@@ -0,0 +1,59 @@
+
+test_that("vfisher gives same results as fisher.test example 1", {
+  TeaTasting <- matrix(
+    c(3, 1, 1, 3), nrow = 2, dimnames = list(
+      Guess = c("Milk", "Tea"),
+      Truth = c("Milk", "Tea")
+  ))
+  fref <- fisher.test(TeaTasting, alternative = "greater")
+  vres <- vfisher.test(
+    TeaTasting[1,1], TeaTasting[1, 2], TeaTasting[2, 1], TeaTasting[2, 2], alternative = "greater"
+  )
+  expect_equal(fref$p.value, vres[, p.value])
+  expect_equal(unname(fref$estimate), vres[, estimate])
+  expect_equal(fref$conf.int[1], vres[, ci.lo])
+  expect_equal(fref$conf.int[2], vres[, ci.hi])
+})
+
+test_that("vfisher gives same results as fisher.test example 2", {
+  Convictions <- matrix(c(2, 10, 15, 3), nrow = 2,
+                        dimnames =
+                          list(c("Dizygotic", "Monozygotic"),
+                               c("Convicted", "Not convicted")))
+  fref <- fisher.test(Convictions, alternative = "less")
+  vres <- vfisher.test(
+    Convictions[1,1], Convictions[1, 2], Convictions[2, 1], Convictions[2, 2], alternative = "less"
+  )
+  expect_equal(fref$p.value, vres[, p.value])
+  expect_equal(unname(fref$estimate), vres[, estimate])
+  expect_equal(fref$conf.int[1], vres[, ci.lo])
+  expect_equal(fref$conf.int[2], vres[, ci.hi])
+})
+
+test_that("vfisher gives same results as fisher.test example 1 & 2", {
+  TeaTasting <- matrix(
+    c(3, 1, 1, 3), nrow = 2, dimnames = list(
+      Guess = c("Milk", "Tea"),
+      Truth = c("Milk", "Tea")
+    ))
+  Convictions <- matrix(c(2, 10, 15, 3), nrow = 2,
+                        dimnames =
+                          list(c("Dizygotic", "Monozygotic"),
+                               c("Convicted", "Not convicted")))
+  fref1 <- fisher.test(TeaTasting)
+  fref2 <- fisher.test(Convictions)
+
+  vres <- vfisher.test(
+    c(TeaTasting[1,1], Convictions[1,1]),
+    c(TeaTasting[1,2], Convictions[1,2]),
+    c(TeaTasting[2,1], Convictions[2,1]),
+    c(TeaTasting[2,2], Convictions[2,2])
+  )
+
+  expect_equal(fref1$p.value, vres[1, p.value])
+  expect_equal(fref2$p.value, vres[2, p.value])
+  expect_equal(unname(fref1$estimate), vres[1, estimate])
+  expect_equal(unname(fref2$estimate), vres[2, estimate])
+  expect_equal(as.numeric(fref1$conf.int), vres[1, c(ci.lo, ci.hi)])
+  expect_equal(as.numeric(fref2$conf.int), vres[2, c(ci.lo, ci.hi)])
+})