# 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

## Sample size search

The sample size returned by `BinarySampleSize()` is checked against a direct scan of the power function.

    check <- function(Test) {
      ss <- BinarySampleSize(0.6, 0.2, 1, 0.025, 0.8, Test)
      data.frame(
        Test = Test,
        N2 = ss$N2,
        power.at.N2 = round(ss$Power, 4),
        power.at.N2.minus.1 = round(
          BinaryPower(0.6, 0.2, ss$N2 - 1, ss$N2 - 1, 0.025, Test)$Power, 4)
      )
    }
    do.call(rbind, lapply(c('Chisq', 'Fisher', 'Fisher-midP', 'Z-pool', 'Boschloo'), check))
    #>          Test N2 power.at.N2 power.at.N2.minus.1
    #> 1       Chisq 23      0.8191              0.7911
    #> 2      Fisher 27      0.8024              0.7993
    #> 3 Fisher-midP 23      0.8155              0.7851
    #> 4      Z-pool 23      0.8088              0.7851
    #> 5    Boschloo 23      0.8088              0.7911

The power reaches the target at the returned sample size and falls short one patient per group below it.

## 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.