Statistical Methods in bbssr

Two-sided conventions

A two-sided version of a discrete conditional test is not unique. Three 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 blaker convention orders the tables by the smaller of their two one-sided tail probabilities and sums the null probabilities of all tables that are at least as extreme as the observed one in this ordering, \[ p(x_{1}, x_{2}) = \sum_{k \,:\, g(k) \le g(x_{1})} f(k), \qquad g(k) = \min\bigl\{\Pr(X_{1} \le k \mid s), \Pr(X_{1} \ge k \mid s)\bigr\} . \] This is formula (2) of Mehrotra, Chan and Berger (2003), and the convention 'blaker' of exact2x2::exact2x2(). The tables of the other tail that it adds form a tail of their own, whose probability is no larger than the tail probability of the observed table. The p-value is therefore at most twice that tail probability and never exceeds the central one. When the two groups are of equal size the conditional distribution is symmetric, and the three conventions give the same p-value.

Fay and Hunsberger (2021, Section 8) illustrate the differences with 8 responders out of 14 against 1 out of 7.

ex <- vapply(c('minlike', 'central', 'blaker'), function(ts) {
  attr(BinaryRR(14, 7, 0.05, 'Fisher', alternative = 'two.sided', tsmethod = ts),
       'p.value')[9, 2]
}, numeric(1))
round(ex, 3)
#> minlike central  blaker 
#>   0.159   0.157   0.087

The table with 4 responders out of 14 against 5 out of 7 is exactly as likely as the observed one, so the minlike convention counts it. Its lower tail probability, 0.0805, exceeds the upper tail probability of the observed table, 0.0783, so the blaker convention leaves it out.

Under the mid-p correction, the minlike and blaker conventions count the tables tied with the observed one in their ordering, the observed one included, with half of their probability, which is the definition of the mid-p value in Fay and Hunsberger (2021, Section 9). The central convention doubles the smaller of the two one-sided mid-p values.

The bound by the central p-value holds for the exact Fisher p-value only. It fails for the mid-p value, and for the Boschloo test, whose ordering of the outcomes changes with the convention, as the two outcomes below show: 5 responders out of 8 against none out of 5 for the mid-p test, and 2 out of 4 against none out of 3 for the Boschloo test.

two_sided_p <- function(N1, N2, Test, x1, x2, ts) {
  attr(BinaryRR(N1, N2, 0.05, Test, alternative = 'two.sided', tsmethod = ts),
       'p.value')[x1 + 1, x2 + 1]
}
data.frame(
  Test = c('Fisher-midP', 'Boschloo'),
  blaker = c(two_sided_p(8, 5, 'Fisher-midP', 5, 0, 'blaker'),
             two_sided_p(4, 3, 'Boschloo', 2, 0, 'blaker')),
  central = c(two_sided_p(8, 5, 'Fisher-midP', 5, 0, 'central'),
              two_sided_p(4, 3, 'Boschloo', 2, 0, 'central'))
)
#>          Test     blaker    central
#> 1 Fisher-midP 0.05361305 0.04351204
#> 2    Boschloo 0.30541895 0.21874043

The central convention has a property that the minlike and blaker conventions lack. 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] TRUE

The minlike and central conventions give different regions of the same nominal size.

minlike <- BinaryRR(N1, N2, 0.05, 'Fisher', alternative = 'two.sided',
                    tsmethod = 'minlike')
central <- BinaryRR(N1, N2, 0.05, 'Fisher', alternative = 'two.sided',
                    tsmethod = 'central')
data.frame(
  tsmethod = c('minlike', 'central'),
  rejected = c(sum(minlike), sum(central)),
  rejected.by.this.only = c(sum(minlike & !central), sum(central & !minlike))
)
#>   tsmethod rejected rejected.by.this.only
#> 1  minlike       22                     2
#> 2  central       20                     0

The chi-squared and Z-pooled tests order outcomes by \(|Z|\) when the alternative is two-sided, so tsmethod does not apply to them.

Refining the maximum between the grid points

The tail probability of each outcome is a polynomial in \(\theta\), and its maximum over the grid is a lower bound of the maximum over the unit interval. A p-value computed on the grid can therefore fall below the exact p-value, and the test can then exceed its nominal level. With ref.pvalue = TRUE the grid is extended by points equally spaced on the arcsine square-root scale, which resolve the narrow local maxima near 0 and 1, and every local maximum on the extended grid is refined by a safeguarded Newton iteration between its two neighbouring grid points. A refined p-value is never smaller than the p-value on the grid.

The size of a fixed-sample design is the largest rejection probability over \(\theta\), which the function below evaluates on a fine grid and refines with optimize() around the largest grid value.

size_fixed <- function(RR) {
  N1 <- attr(RR, 'N1')
  N2 <- attr(RR, 'N2')
  m <- matrix(as.vector(RR), N1 + 1L, N2 + 1L)
  f <- function(t) sum(outer(dbinom(0:N1, N1, t), dbinom(0:N2, N2, t)) * m)
  theta <- seq(0.001, 0.999, by = 0.001)
  v <- vapply(theta, f, numeric(1))
  k <- which.max(v)
  lower <- theta[max(k - 1L, 1L)]
  upper <- theta[min(k + 1L, length(theta))]
  max(v[k], optimize(f, c(lower, upper), maximum = TRUE)$objective)
}
grid.rr <- BinaryRR(32, 32, 0.025, 'Z-pool')
refined.rr <- BinaryRR(32, 32, 0.025, 'Z-pool', ref.pvalue = TRUE)
data.frame(
  p.values = c('grid', 'refined'),
  rejected = c(sum(grid.rr), sum(refined.rr)),
  size = signif(c(size_fixed(grid.rr), size_fixed(refined.rr)), 4)
)
#>   p.values rejected    size
#> 1     grid      342 0.02501
#> 2  refined      340 0.02334

With 32 patients per group the refinement removes 2 outcomes from the rejection region of the Z-pooled test at the one-sided level 0.025, and the size of the test changes from 0.02501 to 0.02334. Every function that computes a rejection region accepts ref.pvalue. The default is FALSE, which keeps the results of earlier versions, and the refinement takes longer than the grid alone.

Ties in the ordering statistic

The 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.05128205

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

p <- attr(BinaryRR(7, 7, 0.025, 'Boschloo'), 'p.value')
c(cell_5_1 = p[6, 2], cell_6_2 = p[7, 3])
#>   cell_5_1   cell_6_2 
#> 0.02867981 0.02867981

The Berger-Boos procedure

Maximizing over the whole unit interval is wasteful, because values of \(\theta\) far from the observed pooled proportion are implausible. Berger and Boos (1994) 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 activates the procedure. Berger and Boos (1994) write \(\beta\) for \(\gamma\) and suggest a small value such as 0.001 or 0.0001. 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.0001087237

Two 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 largest type I error rate

The largest type I error rate is the largest rejection probability on the boundary of the null hypothesis. Without a margin the boundary is \(p_{1} = p_{2} = \theta\). With a margin, on the scale of the risk difference or of the risk ratio, both response probabilities on the boundary are linear functions \(p_{1}(\theta)\) and \(p_{2}(\theta)\) of the pooled response probability \(\theta\), see the bbssr-non-inferiority vignette. BinaryTypeIErrorBSSR() and BinaryAlphaAdjBSSR() find the largest value over an interval of \(\theta\) together with an upper bound, as described in this section.

Let \(x_{11}\) and \(x_{12}\) be the interim numbers of responders, as in the sum above, \(s = x_{11} + x_{12}\) their total, and \(N_{1}(s)\) and \(N_{2}(s)\) the final sample sizes that \(s\) leads to. Given that \(X_{1}\) of the \(N_{1}(s)\) patients of group 1 respond, the number of responders among the first \(n_{11}\) is hypergeometric whatever the response probability, \[ b(x_{11}; n_{11}, p_{1}) \, b\bigl(X_{1} - x_{11}; N_{1}(s) - n_{11}, p_{1}\bigr) = b\bigl(X_{1}; N_{1}(s), p_{1}\bigr) \, h\bigl(x_{11}; n_{11}, N_{1}(s) - n_{11}, X_{1}\bigr) , \] where \(b(x; n, p)\) is the binomial probability and \(h(x; m_{1}, m_{2}, k) = \binom{m_{1}}{x} \binom{m_{2}}{k - x} / \binom{m_{1} + m_{2}}{k}\), and likewise for group 2. Let \(K\) index the distinct pairs \((N_{1K}, N_{2K})\) of final sample sizes and \(\mathcal{R}_{K}\) be the rejection region at these sizes. Summing over the interim outcomes and the final responder counts, the type I error rate is \[ \mathrm{TIE}(\theta) = \sum_{K} \sum_{X_{1} = 0}^{N_{1K}} \sum_{X_{2} = 0}^{N_{2K}} C_{K}(X_{1}, X_{2}) \, b\bigl(X_{1}; N_{1K}, p_{1}(\theta)\bigr) \, b\bigl(X_{2}; N_{2K}, p_{2}(\theta)\bigr) , \] \[ C_{K}(X_{1}, X_{2}) = \mathbf{1}\bigl\{(X_{1}, X_{2}) \in \mathcal{R}_{K}\bigr\} \sum_{(x_{11}, x_{12}) \to K} h(x_{11}; n_{11}, N_{1K} - n_{11}, X_{1}) \, h(x_{12}; n_{12}, N_{2K} - n_{12}, X_{2}) , \] where the inner sum runs over the interim outcomes that lead to \((N_{1K}, N_{2K})\). The weights \(C_{K}\) do not depend on \(\theta\). The fixed-sample design is the case of a single \(K\) without interim patients, with \(C_{1}(X_{1}, X_{2}) = \mathbf{1}\{(X_{1}, X_{2}) \in \mathcal{R}\}\).

On an interval \(\theta_{0} \le \theta \le \theta_{1}\) write \(\theta = (1 - t) \theta_{0} + t \theta_{1}\) with \(0 \le t \le 1\). Then \(p_{j}(\theta) = (1 - t) p_{j}(\theta_{0}) + t p_{j}(\theta_{1})\), and the type I error rate is a polynomial in \(t\), which the package writes in the Bernstein basis \[ B_{k,n}(t) = \binom{n}{k} t^{k} (1 - t)^{n - k} , \qquad k = 0, \ldots, n . \] Three identities give its coefficients without subtraction. First, a response probability \((1 - t) q_{0} + t q_{1}\) is that of a patient who responds with probability \(q_{1}\) with probability \(t\) and with probability \(q_{0}\) otherwise. The number \(J\) of the \(N\) patients of the first kind is binomial with parameters \(N\) and \(t\), so \[ b\bigl(X; N, (1 - t) q_{0} + t q_{1}\bigr) = \sum_{j = 0}^{N} M_{N}(X, j) \, B_{j,N}(t) , \qquad M_{N}(X, j) = \Pr(Y_{j} + Z_{j} = X) , \] where \(Y_{j} \sim \mathrm{Bin}(N - j, q_{0})\) and \(Z_{j} \sim \mathrm{Bin}(j, q_{1})\) are independent, with \(q_{0} = p_{1}(\theta_{0})\) and \(q_{1} = p_{1}(\theta_{1})\) for group 1 and likewise for group 2. Second, the product of two basis polynomials is a multiple of a basis polynomial, \[ B_{j,N_{1}}(t) \, B_{l,N_{2}}(t) = h(j; N_{1}, N_{2}, j + l) \, B_{j + l, N_{1} + N_{2}}(t) , \] so the term of \(K\) is a polynomial of degree \(N_{1K} + N_{2K}\) with the coefficients \[ d_{K,m} = \sum_{j + l = m} h(j; N_{1K}, N_{2K}, m) \sum_{X_{1}, X_{2}} M_{N_{1K}}(X_{1}, j) \, C_{K}(X_{1}, X_{2}) \, M_{N_{2K}}(X_{2}, l) . \] Third, a polynomial of degree \(n\) with the coefficients \(d_{0}, \ldots, d_{n}\) is one of degree \(n + 1\) with the coefficients \[ d'_{i} = \frac{i}{n + 1} \, d_{i - 1} + \Bigl(1 - \frac{i}{n + 1}\Bigr) d_{i} , \qquad i = 0, \ldots, n + 1 , \] where \(d_{-1} = d_{n + 1} = 0\). The terms of the different final sample sizes are raised to the largest degree in this way and added. Every quantity in these steps is non-negative.

The basis polynomials are non-negative and sum to one, so on \(0 \le t \le 1\) the rate lies between the smallest and the largest of its coefficients \(d_{0}, \ldots, d_{n}\), and \(d_{0}\) and \(d_{n}\) are its values at \(t = 0\) and \(t = 1\). The coefficients on the two halves of the interval follow from the algorithm of de Casteljau. With \(d^{(0)}_{i} = d_{i}\) and \(d^{(q)}_{i} = \bigl(d^{(q - 1)}_{i} + d^{(q - 1)}_{i + 1}\bigr) / 2\) for \(i = 0, \ldots, n - q\), the coefficients on the first half are \(d^{(0)}_{0}, d^{(1)}_{0}, \ldots, d^{(n)}_{0}\) and those on the second half are \(d^{(n)}_{0}, d^{(n - 1)}_{1}, \ldots, d^{(0)}_{n}\). The package halves the interval repeatedly and keeps the largest value \(m\) found at the ends of the pieces. A piece is set aside once its largest coefficient is at most \(m + 10^{-12}\), and the largest coefficient \(\bar{m}\) of the pieces set aside is an upper bound of the rate on the whole interval. If fewer than \(10^{5}\) subdivisions suffice, \[ m \le \max_{\theta_{0} \le \theta \le \theta_{1}} \mathrm{TIE}(\theta) \le \bar{m} \le m + 10^{-12} \] up to rounding error. After \(10^{5}\) subdivisions the remaining pieces are set aside with their largest coefficients. The bound \(\bar{m}\) is then still valid, and a warning reports that it can exceed \(m\) by more than \(10^{-12}\).

The interval runs from the smallest to the largest value of theta, restricted with a margin to the values at which both response probabilities lie in the unit interval. BinaryTypeIErrorBSSR() with the default maximize = 'certified' reports \(m\), or the largest value on the grid if that is larger, as TIE, and \(\bar{m}\), or TIE if that is larger by rounding error, as bound. BinaryAlphaAdjBSSR() accepts a level only when \(\bar{m}\) does not exceed the target level. The bbssr-validation vignette compares the result with a dense grid.