Let \(X_{j}\) denote the number of responders in group \(j\), so that \(X_{1} \sim \mathrm{Bin}(N_{1}, p_{1})\) and \(X_{2} \sim \mathrm{Bin}(N_{2}, p_{2})\) independently. Realized counts are written \(x_{1}\) and \(x_{2}\), and \(s = x_{1} + x_{2}\) is the total number of responders. The null hypothesis is \(H_{0}: p_{1} = p_{2}\), and the common value under the null is denoted \(\theta\) and treated as a nuisance parameter.
Every test in the package is defined through its rejection region, a subset of the \((N_{1} + 1) \times (N_{2} + 1)\) grid of possible outcomes. Once that region \(\mathcal{R}\) is available, the power at any pair of response probabilities follows from a single sum,
\[ 1 - \beta = \sum_{(x_{1}, x_{2}) \in \mathcal{R}} \binom{N_{1}}{x_{1}} p_{1}^{x_{1}} (1 - p_{1})^{N_{1} - x_{1}} \binom{N_{2}}{x_{2}} p_{2}^{x_{2}} (1 - p_{2})^{N_{2} - x_{2}} . \]
The type I error rate is the same sum evaluated at \(p_{1} = p_{2} = \theta\), maximized over \(\theta\).
Conditioning on \(s\) removes the nuisance parameter. Under the null the count \(X_{1}\) then follows a hypergeometric distribution, and the one-sided Fisher p-value is
\[ p_{F}(x_{1}, x_{2}) = \Pr(X_{1} \ge x_{1} \mid s) = \sum_{k \ge x_{1}} \frac{\binom{N_{1}}{k} \binom{N_{2}}{s - k}}{\binom{N_{1} + N_{2}}{s}} . \]
The mid-p variant replaces the contribution of the observed table by half of it, giving \(\Pr(X_{1} > x_{1} \mid s) + \tfrac{1}{2} \Pr(X_{1} = x_{1} \mid s)\). This is no longer a valid p-value in the strict sense, so the mid-p test can exceed the nominal level, but it removes much of the conservatism that conditioning introduces.
A two-sided version of a discrete conditional test is not unique. Two
conventions are available through the tsmethod
argument.
The minlike convention sums the null probabilities of
all tables that are no more likely than the observed one, \[
p(x_{1}, x_{2}) = \sum_{k \,:\, f(k) \le f(x_{1})} f(k), \qquad f(k) =
\Pr(X_{1} = k \mid s) .
\] This is the convention of stats::fisher.test.
The central convention doubles the smaller of the two
one-sided tail probabilities and truncates at one, \[
p(x_{1}, x_{2}) = \min\bigl\{1, \; 2 \min(\Pr(X_{1} \le x_{1} \mid s),
\Pr(X_{1} \ge x_{1} \mid s))\bigr\} .
\]
The central convention has a property the minimum-likelihood convention lacks. Its two-sided rejection region at level \(2\alpha\) is exactly the union of the two one-sided rejection regions at level \(\alpha\).
N1 <- 9
N2 <- 7
alpha <- 0.02
two <- BinaryRR(N1, N2, 2 * alpha, 'Fisher',
alternative = 'two.sided', tsmethod = 'central')
upper <- BinaryRR(N1, N2, alpha, 'Fisher')
lower <- t(BinaryRR(N2, N1, alpha, 'Fisher'))
identical(as.vector(two), as.vector(upper | lower))
#> [1] TRUEThe two conventions give different regions of the same nominal size.
data.frame(
tsmethod = c('minlike', 'central'),
rejected = c(
sum(BinaryRR(N1, N2, 0.05, 'Fisher', alternative = 'two.sided',
tsmethod = 'minlike')),
sum(BinaryRR(N1, N2, 0.05, 'Fisher', alternative = 'two.sided',
tsmethod = 'central'))
)
)
#> tsmethod rejected
#> 1 minlike 22
#> 2 central 20The chi-squared and Z-pooled tests order outcomes by \(|Z|\) when the alternative is two-sided, so
tsmethod does not apply to them.
Conditioning is not the only way to eliminate \(\theta\). An exact unconditional test keeps the full binomial model and maximizes the null tail probability over the nuisance parameter,
\[ p(x_{1}, x_{2}) = \sup_{0 \le \theta \le 1} \Pr_{\theta}\bigl(T(X_{1}, X_{2}) \text{ at least as extreme as } T(x_{1}, x_{2})\bigr) , \]
where \(T\) is an ordering statistic. The Z-pooled test uses the two-sample Z statistic with a pooled variance estimator,
\[ Z(x_{1}, x_{2}) = \frac{x_{1} / N_{1} - x_{2} / N_{2}} {\sqrt{\hat{p}(1 - \hat{p})(1 / N_{1} + 1 / N_{2})}}, \qquad \hat{p} = \frac{x_{1} + x_{2}}{N_{1} + N_{2}} , \]
with larger values more extreme. The Boschloo test uses the Fisher p-value itself as the ordering statistic, with smaller values more extreme.
The supremum is approximated by a grid search over \(\theta\). The number of grid points is
controlled by n.grid, which defaults to 100. A finer grid
can only find a larger maximum, so the p-values increase with
n.grid and the test becomes more conservative and more
accurate.
coarse <- attr(BinaryRR(12, 12, 0.025, 'Boschloo', n.grid = 20), 'p.value')
fine <- attr(BinaryRR(12, 12, 0.025, 'Boschloo', n.grid = 2000), 'p.value')
data.frame(all.p.values.increased = all(fine >= coarse - 1e-12),
largest.increase = max(fine - coarse))
#> all.p.values.increased largest.increase
#> 1 TRUE 0.002618907The ordering statistic takes the same value at several outcomes far more often than one might expect. With \(N_{1} = N_{2} = 7\) the outcomes \((x_{1}, x_{2}) = (5, 1)\) and \((6, 2)\) both have a Fisher p-value of \(2/39\).
stat <- attr(BinaryRR(7, 7, 0.025, 'Fisher'), 'p.value')
c(cell_5_1 = stat[6, 2], cell_6_2 = stat[7, 3])
#> cell_5_1 cell_6_2
#> 0.05128205 0.05128205The tail event is defined by “at least as extreme as”, so both outcomes belong to each other’s tail set and must receive the same p-value. Accumulating the null probabilities in an arbitrary order within a tie group would give them different values, and the decision at those outcomes would depend on how the sorting routine happens to break the tie. The package groups tied values explicitly and assigns each group the tail probability accumulated up to its last member.
Maximizing over the whole unit interval is wasteful, because values of \(\theta\) far from the observed pooled proportion are implausible. Berger and Boos proposed maximizing over a \(100(1 - \gamma)\) percent confidence interval \(C_{\gamma}\) for \(\theta\) instead, and adding \(\gamma\) to compensate,
\[ p_{\gamma}(x_{1}, x_{2}) = \gamma + \sup_{\theta \in C_{\gamma}} \Pr_{\theta}\bigl(T(X_{1}, X_{2}) \text{ at least as extreme as } T(x_{1}, x_{2})\bigr) . \]
The package uses an exact Clopper-Pearson interval based on \(s\) responders among \(N_{1} + N_{2}\) patients, so the interval
differs from outcome to outcome. Setting bb.gamma to a
positive value, conventionally \(10^{-4}\) or \(10^{-3}\), activates the procedure. The
resulting test still controls the type I error rate at the nominal
level.
plain <- BinaryRR(15, 15, 0.025, 'Boschloo', n.grid = 200)
bb <- BinaryRR(15, 15, 0.025, 'Boschloo', n.grid = 200, bb.gamma = 1e-4)
p.plain <- attr(plain, 'p.value')
p.bb <- attr(bb, 'p.value')
data.frame(rejected.plain = sum(plain), rejected.berger.boos = sum(bb),
largest.decrease = max(p.plain - p.bb),
largest.increase = max(p.bb - p.plain))
#> rejected.plain rejected.berger.boos largest.decrease largest.increase
#> 1 61 61 0.005134331 0.0001087237Two forces act in opposite directions. Restricting the search lowers the maximum, and the additive \(\gamma\) raises the result. Which one wins varies from outcome to outcome, so the rejection region can grow, shrink or stay as it is. The gain is largest when the outcome is extreme, because the confidence interval for \(\theta\) then excludes the value at which the unrestricted maximum is attained.
The conditional distribution of a p-value satisfies \(\Pr(p_{F} \le c \mid s) \le c\) for any fixed \(c\). Averaging over \(s\) gives \(\Pr_{\theta}(p_{F} \le c) \le c\) for every \(\theta\), so the Boschloo p-value never exceeds the Fisher p-value at the same outcome. The Boschloo rejection region therefore contains the Fisher rejection region, which is the sense in which Boschloo’s test is uniformly more powerful.
fisher <- BinaryRR(15, 15, 0.025, 'Fisher')
boschloo <- BinaryRR(15, 15, 0.025, 'Boschloo', n.grid = 200)
data.frame(rejected.fisher = sum(fisher), rejected.boschloo = sum(boschloo),
fisher.region.contained = all(as.vector(boschloo)[as.vector(fisher)]))
#> rejected.fisher rejected.boschloo fisher.region.contained
#> 1 51 61 TRUEThe type I error rates show how much of the nominal level each test actually spends.
max_type1 <- function(RR, n.grid = 401) {
N1 <- attr(RR, 'N1')
N2 <- attr(RR, 'N2')
m <- matrix(as.vector(RR), N1 + 1L, N2 + 1L)
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)) * m)
}, numeric(1)))
}
tests <- c('Chisq', 'Fisher', 'Fisher-midP', 'Z-pool', 'Boschloo')
t1 <- vapply(tests, function(tst) {
max_type1(BinaryRR(30, 30, 0.025, tst, n.grid = 200), n.grid = 801)
}, numeric(1))
data.frame(
Test = tests, max.type1 = round(t1, 5), exceeds.alpha = t1 > 0.025,
row.names = NULL
)
#> Test max.type1 exceeds.alpha
#> 1 Chisq 0.02770 TRUE
#> 2 Fisher 0.01370 FALSE
#> 3 Fisher-midP 0.02595 TRUE
#> 4 Z-pool 0.02346 FALSE
#> 5 Boschloo 0.02344 FALSEThe Fisher test spends the least, which is the cost of conditioning.
The two unconditional tests spend much more while staying below the
level, which is where their extra power comes from. The chi-squared and
mid-p tests carry no such guarantee, and the exceeds.alpha
column shows what that means at this configuration.
At the interim analysis, \(n_{1}\) and \(n_{2}\) patients have been observed and the total number of responders \(S\) is known. The blinded estimate of the pooled response probability is \(\hat{p} = S / (n_{1} + n_{2})\). With an allocation ratio of \(r\) to 1 and an assumed treatment effect \(\Delta_{A}\), group-specific probabilities are recovered as
\[ \hat{p}_{1} = \hat{p} + \frac{\Delta_{A}}{1 + r}, \qquad \hat{p}_{2} = \hat{p} - \frac{r \Delta_{A}}{1 + r} , \]
truncated to the unit interval. These enter
BinarySampleSize() in place of the original assumptions.
Nothing in this chain requires knowledge of which patient received which
treatment.
Two rules govern what happens next. The unrestricted rule takes the re-estimated sample size as it stands, allowing the trial to end up smaller than planned. The restricted rule raises it to the planned sample size first, so the trial can only grow. The table below applies both to a trial planned at 24 patients per group with an interim analysis at 12.
The unrestricted column is not monotone in \(S\). At a pooled rate of 0.083 the recovered control probability is truncated at zero, which shrinks the recovered risk difference below \(\Delta_{A}\) and pushes the sample size back up.
interim <- data.frame(S = c(2, 4, 6, 8, 10))
interim$pooled <- round(interim$S / 24, 3)
interim$unrestricted <- vapply(interim$S, function(s) {
BinaryBSSR(n1 = 12, n2 = 12, S = s, Delta.A = 0.36, r = 1,
alpha = 0.025, tar.power = 0.8, Test = 'Chisq')$N.final
}, numeric(1))
interim$restricted <- vapply(interim$S, function(s) {
BinaryBSSR(n1 = 12, n2 = 12, S = s, Delta.A = 0.36, r = 1,
alpha = 0.025, tar.power = 0.8, Test = 'Chisq',
restricted = TRUE, N1 = 24, N2 = 24)$N.final
}, numeric(1))
interim
#> S pooled unrestricted restricted
#> 1 2 0.083 40 48
#> 2 4 0.167 30 48
#> 3 6 0.250 40 48
#> 4 8 0.333 50 50
#> 5 10 0.417 58 58BinaryPowerBSSR() evaluates a design by averaging the
conditional power over the distribution of the interim outcome,
\[ 1 - \beta_{\mathrm{BSSR}} = \sum_{x_{1}, x_{2}} \Pr(x_{1} \mid n_{1}, p_{1}) \Pr(x_{2} \mid n_{2}, p_{2}) \, \mathrm{CP}(x_{1}, x_{2}) , \]
where the conditional power \(\mathrm{CP}\) is computed from the rejection region of the final sample size that the interim outcome leads to. The sum runs over every possible interim outcome, so the same rejection region is required many times and is cached.
Boschloo is the default recommendation when the
computation is affordable, since it is exact and uniformly more powerful
than Fisher. Z-pool is close behind and
slightly cheaper. Fisher is exact but conservative, and is
the conventional choice when a regulator expects the classical
procedure. Chisq is useful for exploration and for the
starting value of a sample size search, but does not control the type I
error rate exactly at small sample sizes. Fisher-midP sits
between Fisher and the unconditional tests and is worth
considering when exact control is not a formal requirement.
Berger, R. L. and Boos, D. D. (1994). P values maximized over a confidence set for the nuisance parameter. Journal of the American Statistical Association, 89, 1012-1016.
Boschloo, R. D. (1970). Raised conditional level of significance for the 2x2-table when testing the equality of two probabilities. Statistica Neerlandica, 24, 1-9.
Kieser, M. (2020). Methods and Applications of Sample Size Calculation and Recalculation in Clinical Trials. Springer.
Mehrotra, D. V., Chan, I. S. F. and Berger, R. L. (2003). A cautionary note on exact unconditional inference for a difference between two independent binomial proportions. Biometrics, 59, 441-450.