Validation of bbssr

Gosuke Homma

2026-08-20

library(bbssr)
has_exact <- requireNamespace('Exact', quietly = TRUE)
has_exact2x2 <- requireNamespace('exact2x2', quietly = TRUE)
has_bench <- requireNamespace('microbenchmark', quietly = TRUE)
c(Exact = has_exact, exact2x2 = has_exact2x2, microbenchmark = has_bench)
#>          Exact       exact2x2 microbenchmark 
#>           TRUE           TRUE           TRUE

Chunks that need one of the optional packages are skipped when it is unavailable.

Conditional tests against stats

The Fisher p-values are compared with stats::fisher.test at every outcome of a small grid, for the one-sided alternative and for both two-sided conventions.

N1 <- 7
N2 <- 6
compare_fisher <- function(alternative, tsmethod) {
  got <- attr(BinaryRR(N1, N2, 0.05, 'Fisher', alternative = alternative,
                       tsmethod = tsmethod), 'p.value')
  want <- outer(0:N1, 0:N2, Vectorize(function(i, j) {
    tab <- matrix(c(i, j, N1 - i, N2 - j), nrow = 2)
    if (alternative == 'greater') {
      stats::fisher.test(tab, alternative = 'greater')$p.value
    } else {
      stats::fisher.test(tab)$p.value
    }
  }))
  max(abs(got - want))
}
data.frame(
  comparison = c('one-sided', 'two-sided, minlike'),
  max.absolute.difference = c(
    compare_fisher('greater', 'minlike'),
    compare_fisher('two.sided', 'minlike')
  )
)
#>           comparison max.absolute.difference
#> 1          one-sided            4.440892e-16
#> 2 two-sided, minlike            3.330669e-16

The central convention has no counterpart in stats, so it is checked against its definition.

got <- attr(BinaryRR(N1, N2, 0.05, 'Fisher', alternative = 'two.sided',
                     tsmethod = 'central'), 'p.value')
want <- outer(0:N1, 0:N2, function(i, j) {
  s <- i + j
  pmin(1, 2 * pmin(stats::phyper(i, N1, N2, s),
                   stats::phyper(i - 1, N1, N2, s, lower.tail = FALSE)))
})
max(abs(got - want))
#> [1] 8.881784e-16

Unconditional tests against a direct evaluation

The unconditional p-value is the null tail probability maximized over the nuisance parameter. The reference below forms the tail set of each outcome explicitly and searches the same grid, which is slow but follows the definition with no shared code.

unconditional_ref <- function(stat, N1, N2, n.grid, decreasing) {
  theta <- seq(0, 1, length.out = n.grid)
  joint <- lapply(theta, function(t) outer(dbinom(0:N1, N1, t), dbinom(0:N2, N2, t)))
  p <- matrix(0, nrow = N1 + 1, ncol = N2 + 1)
  for (i in 0:N1) {
    for (j in 0:N2) {
      s0 <- stat[i + 1, j + 1]
      tol <- 1e-10 * pmax(abs(stat), abs(s0))
      mask <- if (decreasing) stat >= s0 - tol else stat <= s0 + tol
      p[i + 1, j + 1] <- max(vapply(joint, function(m) sum(m * mask), numeric(1)))
    }
  }
  pmin(p, 1)
}

N1 <- 8
N2 <- 7
n.grid <- 40
boschloo <- attr(BinaryRR(N1, N2, 0.025, 'Boschloo', n.grid = n.grid), 'p.value')
fisher <- attr(BinaryRR(N1, N2, 0.025, 'Fisher'), 'p.value')
zpool <- attr(BinaryRR(N1, N2, 0.025, 'Z-pool', n.grid = n.grid), 'p.value')
zstats <- outer(0:N1, 0:N2, function(i, j) {
  hat.p <- (i + j) / (N1 + N2)
  z <- (i / N1 - j / N2) / sqrt(hat.p * (1 - hat.p) * (1 / N1 + 1 / N2))
  ifelse(is.finite(z), z, 0)
})
data.frame(
  Test = c('Boschloo', 'Z-pool'),
  max.absolute.difference = c(
    max(abs(boschloo - unconditional_ref(fisher, N1, N2, n.grid, FALSE))),
    max(abs(zpool - unconditional_ref(zstats, N1, N2, n.grid, TRUE)))
  )
)
#>       Test max.absolute.difference
#> 1 Boschloo            2.220446e-16
#> 2   Z-pool            1.110223e-16

Why the tie handling matters

Version 1 of the package accumulated the null probabilities along an arbitrary ordering of the outcomes, which splits groups of tied values. The function below reproduces that behaviour so the two can be compared.

legacy_boschloo <- function(N1, N2, n.grid = 100) {
  stat <- attr(BinaryRR(N1, N2, 0.5, 'Fisher'), 'p.value')
  ord <- order(c(stat))
  x1 <- c(row(stat))[ord] - 1L
  x2 <- c(col(stat))[ord] - 1L
  theta <- seq(0, 1, length.out = n.grid)
  out <- numeric(length(ord))
  for (t in theta) {
    out <- pmax(out, cumsum(dbinom(x1, N1, t) * dbinom(x2, N2, t)))
  }
  p <- stat
  p[ord] <- pmin(1, out)
  p
}

N1 <- 7
N2 <- 7
current <- attr(BinaryRR(N1, N2, 0.025, 'Boschloo'), 'p.value')
legacy <- legacy_boschloo(N1, N2)
c(cells.with.different.p = sum(abs(current - legacy) > 1e-12),
  cells.with.different.decision = sum((current < 0.025) != (legacy < 0.025)))
#>        cells.with.different.p cells.with.different.decision 
#>                            29                             1

The outcomes \((x_{1}, x_{2}) = (5, 1)\) and \((6, 2)\) share a Fisher p-value of \(2/39\), so no test whose ordering statistic is the Fisher p-value can distinguish them.

fisher77 <- attr(BinaryRR(7, 7, 0.5, 'Fisher'), 'p.value')
data.frame(
  outcome = c('x1 = 5, x2 = 1', 'x1 = 6, x2 = 2'),
  fisher.p = c(fisher77[6, 2], fisher77[7, 3]),
  legacy.p = c(legacy[6, 2], legacy[7, 3]),
  current.p = c(current[6, 2], current[7, 3])
)
#>          outcome   fisher.p   legacy.p  current.p
#> 1 x1 = 5, x2 = 1 0.05128205 0.02867981 0.02867981
#> 2 x1 = 6, x2 = 2 0.05128205 0.02161924 0.02867981

The legacy calculation gives the two outcomes different p-values and rejects at one of them but not the other. The rejection region it produces still holds the type I error rate below the nominal level, because the accumulation stops as soon as the running total exceeds alpha, but the decision depends on how the sorting routine breaks the tie rather than on the data.

max_type1 <- function(reject, N1, N2, n.grid = 401) {
  theta <- seq(0, 1, length.out = n.grid)
  max(vapply(theta, function(t) {
    sum(outer(dbinom(0:N1, N1, t), dbinom(0:N2, N2, t)) * reject)
  }, numeric(1)))
}
c(legacy = max_type1(legacy < 0.025, 7, 7),
  current = max_type1(current < 0.025, 7, 7))
#>     legacy    current 
#> 0.02162352 0.01182923

Comparison with the Exact package

N1 <- 10
N2 <- 10
cells <- expand.grid(x1 = c(7, 8, 9), x2 = c(1, 2))
compare_exact <- function(method, Test) {
  ours <- attr(BinaryRR(N1, N2, 0.025, Test, n.grid = 500), 'p.value')
  theirs <- vapply(seq_len(nrow(cells)), function(k) {
    tab <- matrix(c(cells$x1[k], cells$x2[k],
                    N1 - cells$x1[k], N2 - cells$x2[k]), nrow = 2)
    Exact::exact.test(tab, alternative = 'greater', method = method,
                      npNumbers = 500, to.plot = FALSE)$p.value
  }, numeric(1))
  max(abs(ours[cbind(cells$x1 + 1, cells$x2 + 1)] - theirs))
}
data.frame(
  Test = c('Boschloo', 'Z-pool'),
  max.absolute.difference = c(
    compare_exact('boschloo', 'Boschloo'),
    compare_exact('z-pooled', 'Z-pool')
  )
)
#>       Test max.absolute.difference
#> 1 Boschloo            2.109562e-07
#> 2   Z-pool            2.109562e-07

Residual differences come from the grid over the nuisance parameter, which the two packages place differently.

Comparison with exact2x2

N1 <- 10
N2 <- 10
ours <- attr(BinaryRR(N1, N2, 0.05, 'Boschloo', alternative = 'two.sided',
                      tsmethod = 'central', n.grid = 1000), 'p.value')
cells <- list(c(8, 2), c(7, 3), c(9, 1))
data.frame(
  outcome = vapply(cells, function(c) sprintf('x1 = %d, x2 = %d', c[1], c[2]), ''),
  bbssr = vapply(cells, function(c) ours[c[1] + 1, c[2] + 1], numeric(1)),
  exact2x2 = vapply(cells, function(c) {
    exact2x2::boschloo(c[1], N1, c[2], N2, alternative = 'two.sided',
                       tsmethod = 'central')$p.value
  }, numeric(1))
)
#>          outcome       bbssr     exact2x2
#> 1 x1 = 8, x2 = 2 0.011817889 0.0118135656
#> 2 x1 = 7, x2 = 3 0.115318107 0.1152988452
#> 3 x1 = 9, x2 = 1 0.000402448 0.0004021879

Type I error rate of every test

tests <- c('Chisq', 'Fisher', 'Fisher-midP', 'Z-pool', 'Boschloo')
type1 <- function(Test, alternative, N = 30, alpha = 0.025) {
  RR <- BinaryRR(N, N, alpha, Test, alternative = alternative, n.grid = 200)
  max_type1(matrix(as.vector(RR), N + 1, N + 1), N, N, n.grid = 801)
}
one <- vapply(tests, type1, numeric(1), 'greater')
two <- vapply(tests, type1, numeric(1), 'two.sided')
data.frame(
  Test = tests,
  one.sided = round(one, 5),
  two.sided = round(two, 5),
  exceeds.alpha = one > 0.025 | two > 0.025,
  row.names = NULL
)
#>          Test one.sided two.sided exceeds.alpha
#> 1       Chisq   0.02770   0.02741          TRUE
#> 2      Fisher   0.01370   0.01350         FALSE
#> 3 Fisher-midP   0.02595   0.02735          TRUE
#> 4      Z-pool   0.02346   0.02293         FALSE
#> 5    Boschloo   0.02344   0.02267         FALSE

The nominal level is 0.025. The three exact tests stay below it for both alternatives, at every sample size and every value of the nuisance parameter. The chi-squared and mid-p tests carry no such guarantee, and the exceeds.alpha column records where they overshoot at this configuration.

Speed

The whole rejection region is computed in one pass, so obtaining a power curve costs little more than a single p-value. The inner loop over the nuisance parameter runs in compiled code.

microbenchmark::microbenchmark(
  Chisq = BinaryRR(50, 50, 0.025, 'Chisq'),
  Fisher = BinaryRR(50, 50, 0.025, 'Fisher'),
  `Z-pool` = BinaryRR(50, 50, 0.025, 'Z-pool'),
  Boschloo = BinaryRR(50, 50, 0.025, 'Boschloo'),
  `Boschloo, Berger-Boos` = BinaryRR(50, 50, 0.025, 'Boschloo', bb.gamma = 1e-4),
  times = 10L,
  unit = 'ms'
)
#> Unit: milliseconds
#>                   expr     min      lq     mean   median      uq     max neval
#>                  Chisq  0.6408  1.0607  1.13466  1.11495  1.3027  1.5610    10
#>                 Fisher  3.4122  3.9542  4.89522  4.83035  5.5202  7.0450    10
#>                 Z-pool  5.3481  7.4868  7.47379  7.73120  7.8269  8.1077    10
#>               Boschloo  8.0673 11.1210 12.15003 12.71815 13.4889 15.5809    10
#>  Boschloo, Berger-Boos 20.3736 26.6152 26.50769 27.05400 27.6324 28.3149    10
#>    cld
#>  a    
#>   b   
#>    c  
#>     d 
#>      e

The Berger-Boos variant costs roughly twice as much as the plain Boschloo test, because the confidence bounds of every possible responder total are added to the grid over the nuisance parameter.

A cell-by-cell comparison against a package that computes one p-value at a time is not like for like, since BinaryRR() returns all 2601 p-values of the grid above. The comparison below therefore charges Exact only for the p-values of a single row.

N1 <- 20
N2 <- 20
microbenchmark::microbenchmark(
  bbssr.whole.grid = BinaryRR(N1, N2, 0.025, 'Boschloo'),
  Exact.one.row = for (x1 in 0:N1) {
    Exact::exact.test(matrix(c(x1, 5, N1 - x1, N2 - 5), nrow = 2),
                      alternative = 'greater', method = 'boschloo',
                      npNumbers = 100, to.plot = FALSE)
  },
  times = 5L,
  unit = 'ms'
)
#> Unit: milliseconds
#>              expr      min       lq      mean   median       uq      max neval
#>  bbssr.whole.grid   2.3065   3.4072   3.42304   3.5341   3.6729   4.1945     5
#>     Exact.one.row 300.8874 317.9722 333.91644 322.4450 360.3704 367.9072     5
#>  cld
#>   a 
#>    b

Summary

The conditional tests reproduce stats::fisher.test to machine precision. The unconditional tests reproduce a direct evaluation of their definition to machine precision, and agree with Exact and exact2x2 up to the difference between the two grids over the nuisance parameter. The three exact tests hold the type I error rate below the nominal level for one-sided and two-sided alternatives alike.