--- title: "Statistical Methods in bbssr" author: "Gosuke Homma" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{Statistical Methods in bbssr} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, fig.align = "center" ) ``` ```{r load} library(bbssr) ``` ## Notation 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$. ## Conditional tests 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. ### Two-sided conventions 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$. ```{r central-identity} 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)) ``` The two conventions give different regions of the same nominal size. ```{r ts-compare} 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')) ) ) ``` The chi-squared and Z-pooled tests order outcomes by $|Z|$ when the alternative is two-sided, so `tsmethod` does not apply to them. ## Unconditional tests 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. ```{r grid} 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)) ``` ### 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$. ```{r ties} stat <- attr(BinaryRR(7, 7, 0.025, 'Fisher'), 'p.value') c(cell_5_1 = stat[6, 2], cell_6_2 = stat[7, 3]) ``` 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. ```{r ties-equal} p <- attr(BinaryRR(7, 7, 0.025, 'Boschloo'), 'p.value') c(cell_5_1 = p[6, 2], cell_6_2 = p[7, 3]) ``` ### The Berger-Boos refinement 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. ```{r berger-boos} 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)) ``` 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. ## Relationships between the tests 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. ```{r nesting} 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)])) ``` The type I error rates show how much of the nominal level each test actually spends. ```{r type1} 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 ) ``` The 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. ## Blinded sample size re-estimation 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. ```{r bssr-rules} 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 ``` `BinaryPowerBSSR()` 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. ## Choosing a test `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. ## References 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.