--- title: "Type I Error Rate of a Re-estimation Design" author: "Gosuke Homma" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{Type I Error Rate of a Re-estimation Design} %\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) ``` ## The type I error rate of a re-estimation design With blinded sample size re-estimation the final sample size depends on the pooled number of interim responders, whose distribution depends on the common response probability $\theta$. The type I error rate of a re-estimation design is computed for the design itself, by summing the rejection probability over the interim outcomes and the outcomes of the second stage. Under the null hypothesis $p_{1} = p_{2} = \theta$. Let $n_{11}$ and $n_{12}$ be the interim sample sizes of the two groups, $s$ the pooled number of interim responders, $n_{21}(s)$ and $n_{22}(s)$ the second-stage sample sizes reached from $s$, and $\mathcal{R}_{s}$ the rejection region at the final sample size reached from $s$. The type I error rate is $$ \mathrm{TIE}(\theta) = \sum_{s = 0}^{n_{11} + n_{12}} \Pr_{\theta}(S = s) \, \Pr_{\theta}\bigl((X_{1}, X_{2}) \in \mathcal{R}_{s} \mid S = s\bigr) , $$ a polynomial in $\theta$ that the package evaluates exactly. ## The rate over the nuisance parameter `BinaryTypeIErrorBSSR()` evaluates the type I error rate of a re-estimation design and of the fixed-sample design with the initial sample size over a grid of $\theta$, by default `seq(0, 1, by = 0.005)`. The largest value is not read off the grid. With the default `maximize = 'certified'` the function finds the largest value over the whole interval spanned by the grid, here $[0, 1]$, together with an upper bound that the rate does not exceed anywhere in that interval. The `bbssr-statistical-methods` vignette describes the method. The design below has an assumed difference of 0.3, an initial sample size of 39 per group, an interim analysis after 20 per group, re-estimation by formula (21.3) of Kieser (2020) and the chi-squared test. ```{r tie} design <- list(Delta.A = 0.3, N1 = 39, N2 = 39, n.interim = c(20, 20), r = 1, alpha = 0.025, tar.power = 0.8, Test = 'Chisq', ss.method = 'standard') tie <- do.call(BinaryTypeIErrorBSSR, design) tie attr(tie, 'max') plot(tie) ``` ```{r tie-max, include = FALSE} m <- attr(tie, 'max') stopifnot(nrow(m) == 2, all(m$bound >= m$TIE), all(m$bound - m$TIE <= 1e-12 + 1e-15)) ``` The largest type I error rate is `r signif(m$TIE[1], 4)` at $\theta =$ `r signif(m$theta[1], 4)` for the re-estimation design and `r signif(m$TIE[2], 4)` at $\theta =$ `r signif(m$theta[2], 4)` for the fixed-sample design. The column `bound` of the attribute `max` exceeds these values by at most $10^{-12}$, so the largest rates are known to that accuracy whatever the spacing of the grid. The chi-squared test is not exact, so the fixed-sample design can already exceed the nominal level, and the re-estimation changes the rate again. ## Exact tests in the final analysis The same design with an exact test in the final analysis gives the following largest rates. ```{r exact-tests} exact.max <- do.call(rbind, lapply(c('Fisher', 'Boschloo'), function(tst) { m <- attr(do.call(BinaryTypeIErrorBSSR, modifyList(design, list(Test = tst))), 'max') data.frame(Test = tst, Design = m$Design, theta = signif(m$theta, 4), TIE = signif(m$TIE, 4)) })) exact.max ``` ```{r exact-tests-check, include = FALSE} stopifnot(all(exact.max$TIE[exact.max$Design == 'BSSR'] <= 0.025)) ``` For this rule the largest rate of the re-estimation design is `r exact.max$TIE[exact.max$Test == 'Fisher' & exact.max$Design == 'BSSR']` with Fisher's exact test and `r exact.max$TIE[exact.max$Test == 'Boschloo' & exact.max$Design == 'BSSR']` with the Boschloo test, against the nominal 0.025. ## The adjusted significance level When the rate exceeds the nominal level, `BinaryAlphaAdjBSSR()` lowers the nominal level until the largest type I error rate does not exceed the target level, following Kieser and Friede (2000) and Friede and Kieser (2004). With `adjust = 'test'` (the default) the adjusted level is applied to the final analysis only. The re-estimated sample sizes then do not depend on the level, the largest rate is a non-decreasing step function of it, and the level is found by bisection. Every level between the largest p-value that is rejected and the smallest p-value that is not gives the same rejection regions, and the reported level is the largest value with at most six significant digits in that range, so it can be used as printed. Both searches below use the grid `theta.grid` from 0.01 to 0.99 in steps of 0.01, which is coarser than the default to keep them short. The grid serves the search: a level is assessed on the grid, and with the default `maximize = 'certified'` the level found is accepted only if the upper bound of the type I error rate over the whole interval $[0.01, 0.99]$ does not exceed the target level. The columns `max.TIE.bound` and `max.TIE.adj.bound` of the result hold the bounds at the nominal and at the adjusted level. ```{r adjust-test} theta.grid <- seq(0.01, 0.99, by = 0.01) adj <- do.call(BinaryAlphaAdjBSSR, c(design, list(theta = theta.grid))) adj ``` With `adjust = 'both'` the adjusted level is used in the re-estimation as well. The re-estimated sample sizes then change with the level, the rate need not be monotone in it, and the level is lowered in steps of `step` until the rate is controlled. The result is the first level on this grid of steps at which the rate is controlled. The levels that control the rate need not form an interval, so the result can change with the step. The published values of Friede and Kieser (2004, Section 5) and Kieser (2020, Example 21.1) are reproduced only when the re-estimation also uses the adjusted level, as the script `reproduce-published.R` shows (see the `bbssr-validation` vignette). The search below uses a step of 0.00005, coarser than the default, to keep it short. ```{r adjust-both} adj.both <- do.call(BinaryAlphaAdjBSSR, c(design, list(adjust = 'both', step = 5e-5, theta = theta.grid))) adj.both ``` The adjusted level of the re-estimation design is `r adj$alpha.adj[1]` when it applies to the final analysis only and `r adj.both$alpha.adj[1]` when it also applies to the re-estimation. At these levels the upper bounds of the type I error rate over $[0.01, 0.99]$ are `r signif(adj$max.TIE.adj.bound[1], 4)` and `r signif(adj.both$max.TIE.adj.bound[1], 4)`. The power at the adjusted level follows from `BinaryPowerBSSR()` with `alpha` set to the adjusted level. The re-estimation keeps the target level through `ss.alpha` when the level was adjusted with `adjust = 'test'`. ```{r adjusted-power} adjusted <- modifyList(design, list(alpha = adj$alpha.adj[1], ss.alpha = 0.025)) power.adj <- do.call(BinaryPowerBSSR, c(adjusted, list(p = c(0.3, 0.4, 0.5), Delta.T = 0.3))) power.nominal <- do.call(BinaryPowerBSSR, c(design, list(p = c(0.3, 0.4, 0.5), Delta.T = 0.3))) data.frame(p = power.adj$p, power.nominal = round(power.nominal$power.BSSR, 4), power.adjusted = round(power.adj$power.BSSR, 4)) ``` ## Control within a confidence interval The largest type I error rate over the whole range of the nuisance parameter can occur at values that the data of the completed trial make implausible. Kieser (2020, Section 22.1.3) describes another way to control the type I error rate. A value $0 < \gamma < \alpha$ is fixed in advance, the $(1 - \gamma)$ confidence interval for the nuisance parameter is computed from the data of the completed trial, and the final analysis uses a level at which the largest type I error rate over this interval does not exceed $\alpha - \gamma$. The type I error rate is then controlled at $\alpha$. `BinaryAlphaAdjBSSR()` finds this level when `theta` spans the interval and `alpha` is set to $\alpha - \gamma$. The interval is known only at the end of the trial, so the re-estimation keeps the nominal level, through `adjust = 'test'` and `ss.alpha`. Example 23.1 of Kieser (2020) applies the approach to the FreezeAF trial, a non-inferiority trial with the margin 0.15 and the test of Farrington and Manning (1990), equal assumed rates of 0.78, an initial sample size of 122 per group and a pilot of 100 patients. The rate is evaluated on the boundary of the null hypothesis against the pooled response probability, as described in the `bbssr-non-inferiority` vignette, on the grid of the book with steps of 0.005. For control over the whole range, the book applies the adjusted level to the re-estimation as well, which is `adjust = 'both'`. A step of 0.0001 keeps the search short. ```{r ci-whole} freeze <- list(Delta.A = 0, N1 = 122, N2 = 122, n.interim = c(50, 50), r = 1, alpha = 0.025, tar.power = 0.8, Test = 'Farrington-Manning', ss.method = 'standard', margin = 0.15) adj.whole <- do.call(BinaryAlphaAdjBSSR, c(freeze, list(theta = seq(0.075, 0.925, by = 0.005), maximize = 'grid', adjust = 'both', step = 1e-4))) adj.whole[1, c('Design', 'max.TIE', 'alpha.adj', 'max.TIE.adj', 'theta.adj')] ``` The trial observed 210 responders among 291 patients. With $\gamma = 0.0001$ the interval is the 99.99 per cent Clopper-Pearson interval, which `binom.test()` returns. The book calls it a 95 per cent interval, but the interval it gives is the 99.99 per cent interval. The code below computes the largest type I error rate within the interval at the nominal level and the adjusted level. ```{r ci-interval} gamma <- 1e-4 ci <- binom.test(210, 291, conf.level = 1 - gamma)$conf.int ci.grid <- seq(ci[1], ci[2], length.out = 42) tie.ci <- do.call(BinaryTypeIErrorBSSR, c(freeze, list(theta = ci.grid, maximize = 'grid'))) adj.ci <- do.call(BinaryAlphaAdjBSSR, c(modifyList(freeze, list(alpha = 0.025 - gamma)), list(ss.alpha = 0.025, theta = ci.grid, maximize = 'grid'))) adj.ci[1, c('Design', 'alpha', 'alpha.adj', 'max.TIE.adj', 'theta.adj')] ``` Within the interval $[`r round(ci[1], 3)`, `r round(ci[2], 3)`]$ the largest rate of the re-estimation design is `r signif(max(tie.ci$TIE.BSSR), 3)`, only slightly below the largest rate `r signif(adj.whole$max.TIE[1], 3)` over the whole range, as the book notes. The adjusted level that keeps the rate within the interval at or below $\alpha - \gamma = `r 0.025 - gamma`$ is `r signif(adj.ci$alpha.adj[1], 3)`, and the level that controls the rate over the whole range is `r adj.whole$alpha.adj[1]`. The book reports 0.0241 for both, and `reproduce-published.R` reproduces these values with a step of $10^{-5}$ (see the `bbssr-validation` vignette). ## Conditional rejection probabilities `BinaryCondRejectBSSR()` takes the calculation apart. Given the pooled number $s$ of interim responders and the pooled number $s_{2}$ of second-stage responders, the number of responders of group 1 in each stage is hypergeometric under the null hypothesis, and the two stages are independent. The rejection probability given $s$ and $s_{2}$, $$ \mathrm{CRP}(s, s_{2}) = \sum_{x_{11}} \sum_{x_{21}} h_{1}(x_{11} \mid s) \, h_{2}(x_{21} \mid s_{2}) \, \mathbf{1}\bigl\{(x_{11} + x_{21}, \, s + s_{2} - x_{11} - x_{21}) \in \mathcal{R}_{s}\bigr\} , $$ where $h_{1}$ and $h_{2}$ are the hypergeometric probabilities of the two stages, does not depend on $\theta$. The type I error rate is its average over the binomial distributions of $s$ and $s_{2}$, $$ \mathrm{TIE}(\theta) = \sum_{s} \sum_{s_{2}} b(s; n_{11} + n_{12}, \theta) \, b\bigl(s_{2}; n_{21}(s) + n_{22}(s), \theta\bigr) \, \mathrm{CRP}(s, s_{2}) , $$ where $b$ is the binomial probability. The column `CRP.total` is the rejection probability given only the total $s + s_{2}$, with the responder count of group 1 treated as one hypergeometric count among all final patients. ```{r crp} crp <- BinaryCondRejectBSSR( Delta.A = 0.3, N1 = 12, N2 = 12, n.interim = c(6, 6), r = 1, alpha = 0.025, tar.power = 0.8, Test = 'Chisq', ss.method = 'standard', theta = c(0.3, 0.5) ) crp plot(crp) ``` The decomposition by $s$ is returned when `theta` is given. ```{r by-s} by.s <- attr(crp, 'by.s') head(by.s) attr(crp, 'TIE') plot(crp, what = 'by.s') ``` `TIE.s` is the type I error rate given $s$, and `contribution` is its product with the probability `prob.s` of $s$. The contributions add up to the type I error rate at each value of `theta`, which is `r paste(signif(attr(crp, 'TIE')$TIE, 4), collapse = ' and ')` here. ## Non-inferiority With a non-inferiority margin, `BinaryTypeIErrorBSSR()` and `BinaryAlphaAdjBSSR()` evaluate the rate on the boundary of the null hypothesis, parametrized by the pooled response probability, see the `bbssr-non-inferiority` vignette. `BinaryCondRejectBSSR()` covers tests of superiority only, since on the boundary of a non-inferiority hypothesis the two response probabilities differ and the conditional distributions of the responder counts depend on them. ## References Farrington, C. P. and Manning, G. (1990). Test statistics and sample size formulae for comparative binomial trials with null hypothesis of non-zero risk difference or non-unity relative risk. *Statistics in Medicine*, 9, 1447-1454. Friede, T. and Kieser, M. (2004). Sample size recalculation for binary data in internal pilot study designs. *Pharmaceutical Statistics*, 3, 269-279. Kieser, M. (2020). *Methods and Applications of Sample Size Calculation and Recalculation in Clinical Trials*. Springer. Kieser, M. and Friede, T. (2000). Re-calculating the sample size in internal pilot study designs with control of the type I error rate. *Statistics in Medicine*, 19, 901-911.