---
title: "Introduction to bbssr"
author: "Gosuke Homma"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Introduction to 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)
```

## What the package does

A trial with a binary endpoint needs an assumed response probability in each group before the first patient is enrolled. If the pooled response probability turns out to differ from what was assumed, the trial ends up either underpowered or larger than it needed to be. Blinded sample size re-estimation addresses this by looking at the pooled number of responders partway through the trial and adjusting the remaining enrolment, without ever splitting the data by treatment group.

The package covers the whole chain of calculations that this requires. Everything rests on `BinaryRR()`, which returns the rejection region of an exact test over the grid of possible outcomes. `BinaryPower()` sums the binomial probability mass over that region, `BinarySampleSize()` searches for the smallest sample size attaining a target power, `BinaryPowerBSSR()` averages the power over the distribution of the interim data, and `BinaryBSSR()` applies the re-estimation to a single observed interim data set.

## The five tests

```{r tests}
tests <- c('Chisq', 'Fisher', 'Fisher-midP', 'Z-pool', 'Boschloo')
```

`Chisq` is the Pearson chi-squared test without a continuity correction. `Fisher` is the Fisher exact test and `Fisher-midP` its mid-p variant. `Z-pool` and `Boschloo` are exact unconditional tests, which remove the nuisance parameter by maximizing the null tail probability over the common response probability rather than by conditioning on the observed margin.

The unconditional tests are the most powerful of the five and the slowest to compute. Only `Fisher`, `Z-pool` and `Boschloo` guarantee that the type I error rate stays below the nominal level for every value of the nuisance parameter. The chi-squared and mid-p tests trade that guarantee for shorter rejection thresholds.

## Rejection regions

```{r rr}
RR <- BinaryRR(N1 = 10, N2 = 10, alpha = 0.025, Test = 'Boschloo')
RR
```

```{r rr-plot}
plot(RR)
```

The region is monotone: more responders in group 1, or fewer in group 2, can only move an outcome into the rejection region. Because the Boschloo p-value never exceeds the Fisher p-value, the Boschloo region always contains the Fisher region.

```{r rr-compare}
n_rejected <- vapply(tests, function(tst) {
  sum(BinaryRR(10, 10, 0.025, tst))
}, numeric(1))
n_rejected
```

## Power and sample size

```{r power}
BinaryPower(p1 = 0.6, p2 = 0.3, N1 = 40, N2 = 40, alpha = 0.025, Test = 'Fisher')
```

```{r power-curve}
pw <- BinaryPower(p1 = seq(0.35, 0.75, by = 0.05), p2 = rep(0.3, 9),
                  N1 = 40, N2 = 40, alpha = 0.025, Test = 'Fisher')
plot(pw)
```

```{r samplesize}
ss <- BinarySampleSize(p1 = 0.6, p2 = 0.3, r = 1, alpha = 0.025,
                       tar.power = 0.8, Test = 'Fisher')
ss
```

The exact power is not monotone in the sample size. Adding one patient changes the set of attainable significance levels, and the largest level below alpha can fall rather than rise. The following table shows the effect around the selected sample size.

```{r sawtooth}
n2 <- (ss$N2 - 6):(ss$N2 + 5)
data.frame(
  N2 = n2,
  Power = round(vapply(n2, function(n) {
    BinaryPower(0.6, 0.3, n, n, 0.025, 'Fisher')$Power
  }, numeric(1)), 4)
)
```

This is why the sample size search evaluates the exact power at every candidate rather than inverting a smooth approximation.

## Blinded sample size re-estimation

Consider a trial designed for a response probability of 0.45 in group 1 and 0.09 in group 2, so an assumed treatment effect of 0.36 and a pooled probability of 0.27.

The variance of a binary endpoint is largest at one half, so at a fixed risk difference the required sample size grows as the pooled probability approaches that value. A trial sized at a pooled probability of 0.27 is therefore over-powered if the true value turns out lower, and under-powered if it turns out higher.

```{r bssr-plan}
BinarySampleSize(p1 = 0.45, p2 = 0.09, r = 1, alpha = 0.025,
                 tar.power = 0.8, Test = 'Z-pool')
```

`BinaryPowerBSSR()` traces both effects. The scenarios below run from a pooled probability of 0.19 up to 0.37, with the design fixed at the sample size above.

```{r bssr-power}
res <- BinaryPowerBSSR(
  p = seq(0.19, 0.37, by = 0.03),
  Delta.A = 0.36, Delta.T = 0.36,
  N1 = 24, N2 = 24, omega = 0.5, r = 1,
  alpha = 0.025, tar.power = 0.8, Test = 'Z-pool'
)
res
```

```{r bssr-plot}
plot(res)
```

The `power.TRAD` column falls steadily as the pooled probability rises, from well above the target at the low end to well below it at the high end. The `power.BSSR` column stays nearer the target at both ends, since the re-estimation removes patients in the first case and adds them in the second. The `E.N` column shows the expected total sample size that the re-estimation settles on, against the 48 patients of the fixed design.

## Restricted and unrestricted rules

Under the unrestricted rule the trial may end up smaller than planned when the interim data suggest that fewer patients suffice. Under the restricted rule the planned sample size acts as a floor, so the trial can only grow.

```{r rules}
args <- list(
  p = seq(0.19, 0.37, by = 0.06),
  Delta.A = 0.36, Delta.T = 0.36, N1 = 24, N2 = 24, omega = 0.5, r = 1,
  alpha = 0.025, tar.power = 0.8, Test = 'Chisq'
)
unrestricted <- do.call(BinaryPowerBSSR, c(args, list(restricted = FALSE)))
restricted <- do.call(BinaryPowerBSSR, c(args, list(restricted = TRUE)))
data.frame(
  p = unrestricted$p,
  power.unrestricted = round(unrestricted$power.BSSR, 4),
  power.restricted = round(restricted$power.BSSR, 4),
  EN.unrestricted = round(unrestricted$E.N, 1),
  EN.restricted = round(restricted$E.N, 1)
)
```

For any single interim outcome that calls for more patients than planned the two rules give the same answer, since neither caps the increase. They part company when the interim outcome calls for fewer. The columns below are expectations over all interim outcomes, so they converge rather than coincide as the pooled probability rises. The unrestricted rule then takes the saving, ending with a smaller trial and a power closer to the target. The restricted rule keeps the planned sample size, which leaves the trial over-powered but never smaller than what the protocol promised.

## Re-estimating from observed interim data

The functions above describe a design before the trial starts. Once the trial is running, the interim analysis produces one number that can be shared without unblinding, namely the total count of responders. `BinaryBSSR()` turns that number into an enrolment decision.

```{r bssr-interim}
BinaryBSSR(n1 = 12, n2 = 12, S = 8, Delta.A = 0.36, r = 1,
           alpha = 0.025, tar.power = 0.8, Test = 'Z-pool')
```

The interim analysis sits at half of the planned 24 patients per group. Eight responders among 24 patients give a blinded pooled rate of 0.33, above the assumed 0.27, so the re-estimated total exceeds the plan and the second stage has to enrol more patients than originally scheduled. The `bbssr-interim-reestimation` vignette works through this in more detail.

## Two-sided tests

Every test accepts `alternative = 'two.sided'`. For the conditional tests the two-sided p-value follows one of two conventions, selected with `tsmethod`.

```{r two-sided}
pw <- function(alpha, alternative, tsmethod = 'minlike') {
  BinaryPower(0.6, 0.3, 40, 40, alpha, 'Fisher',
              alternative = alternative, tsmethod = tsmethod)$Power
}
round(c(
  one.sided.alpha.0.025 = pw(0.025, 'greater'),
  two.sided.alpha.0.025.minlike = pw(0.025, 'two.sided', 'minlike'),
  two.sided.alpha.0.025.central = pw(0.025, 'two.sided', 'central'),
  two.sided.alpha.0.05.minlike = pw(0.05, 'two.sided', 'minlike'),
  two.sided.alpha.0.05.central = pw(0.05, 'two.sided', 'central')
), 4)
```

At the same nominal level the two-sided test is the less powerful, since half of the level is spent on a direction the alternative does not point in. At twice the level it comes back to the one-sided test, because the lower tail contributes almost nothing to the power when the true effect is positive.

The `bbssr-statistical-methods` vignette explains the two conventions and the Berger-Boos refinement of the unconditional tests.

## Where to go next

The `bbssr-statistical-methods` vignette derives the tests and the re-estimation rules. The `bbssr-interim-reestimation` vignette follows a single trial from planning to the interim decision. The `bbssr-validation` vignette compares the package against `stats`, `Exact` and `exact2x2`, and reports timings.