#' ---
#' title: "Defining group-sequential boundaries with rpact"
#' author: "Marcel Wolbers, Gernot Wassmer, and Friedrich Pahlke"
#' date: "Last change: `r format(Sys.time(), '%d %B, %Y')`"
#' number: 1
#' header-includes:
#' - \usepackage{fancyhdr}
#' - \pagestyle{fancy}
#' - \setlength{\headheight}{23pt}
#' - \fancyfoot[C]{www.rpact.com}
#' - \fancyfoot[L]{\thepage}
#' output:
#' rmarkdown::html_document:
#' highlight: pygments
#' number_sections: yes
#' self_contained: yes
#' toc: yes
#' toc_depth: 3
#' toc_float: yes
#' css: style.css
#' includes:
#' before_body: header.html
#' after_body: footer2.html
#' ---
#'
#'
#'
#' # Summary {-}
#'
#' This [R Markdown](https://rmarkdown.rstudio.com) document provides
#' example code for the the definition of the most commonly used group-sequential
#' boundaries in [rpact](https://cran.r-project.org/package=rpact).
#'
#' # Introduction
#'
#' In rpact, **sample size calculation for a group-sequential trial proceeds by
#' following the same two steps regardless whether the endpoint is a continuous,
#' binary, or a time-to-event endpoint**:
#'
#' 1. **Define the (abstract) group-sequential boundaries** of the design using
#' the function `getDesignGroupSequential()`.
#' 2. **Calculate sample size for the endpoint of interest** by feeding the
#' abstract boundaries from step 1. into specific functions for the endpoint of
#' interest. This step uses functions such as `getSampleSizeMeans()` (for
#' continuous endpoints), `getSampleSizeRates()` (for binary endpoints), and
#' `getSampleSizeSurvival()` (for survival endpoints).
#'
#' The mathematical rationale for this two-step approach is that all
#' group-sequential trials, regardless of the chosen endpoint type, rely on the
#' fact that the $z$-scores at different interim stages follow the same "canonical
#' joint multivariate distribution" (at least asymptotically).
#'
#' This document covers the more abstract first step,
#' **Step 2 is not covered in this document but it is covered in the separate
#' endpoint-specific R Markdown files for continuous, binary, and time to event
#' endpoints.** Of note, step 1 can be omitted for trials without interim
#' analyses.
#'
#' These examples are not intended to replace the official rpact documentation and
#' help pages but rather to supplement them.
#'
#' In general, rpact supports both one-sided and two-sided group-sequential
#' designs. If futility boundaries are specified, however, only one-sided tests are
#' permitted. **For simplicity, it is often preferred to use one-sided tests for
#' group-sequential designs** (typically, with $\alpha = 0.025$).
#'
#' **First, load the rpact package**
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
library(rpact)
packageVersion("rpact") # version should be version 2.0.5 or later
#'
#'
#' # Designs with efficacy interim analyses
#'
#' ## O'Brien & Fleming type $\alpha$-spending
#'
#' **Example:**
#'
#' - Interim analyses at information fractions 33%, 67%, and 100% (
#' `informationRates = c(0.33, 0.67, 1)` ).
#' [Note: For equally spaced interim analyses, one can also specify
#' the maximum number of stages (`kMax`, including the final analysis)
#' instead of the `informationRates`.]
#' - Lan & DeMets $\alpha$-spending approximation to the O'Brien & Fleming boundary
#' (`typeOfDesign = "asOF"`)
#' - $\alpha$-spending approaches allow for flexible timing of interim analyses and
#' corresponding adjustment of boundaries.
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
design <- getDesignGroupSequential(sided = 1, alpha = 0.025,
informationRates = c(0.33, 0.67, 1), typeOfDesign = "asOF")
#'
#' ## Standard O'Brien & Fleming boundary
#'
#' The originally published O'Brien & Fleming boundary is obtained via `typeOfDesign = "OF"` which
#' is also the default (therefore, if you do not specify `typeOfDesign`, this boundary is selected).
#' Note that strict Type I error control is only guaranteed for standard boundaries
#' without $\alpha$-spending if the pre-defined interim schedule (i.e., the
#' information fractions at which interim analyses are conducted) is exactly adhered to.
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
design <- getDesignGroupSequential(sided = 1, alpha = 0.025,
informationRates = c(0.33, 0.67, 1), typeOfDesign = "OF")
#'
#'
#' ## Standard Pocock and Haybittle & Peto boundaries
#' Pocock (`typeOfDesign = "P"` for standard boundary, `typeOfDesig = "asP"` for
#' $\alpha$-spending version) or Haybittle & Peto (`typeOfDesign = "HP"`) boundaries
#' (reject at interim if $z$-value exceeds 3) is obtained with, for example,
## ---- include = TRUE, echo = TRUE---------------------------------------------
design <- getDesignGroupSequential(sided = 1, alpha = 0.025,
informationRates = c(0.33,0.67,1), typeOfDesign = "P")
#'
#' ## Other pre-defined boundaries
#' - Kim & DeMets $\alpha$-spending (`typeOfDesign = "asKD`) with parameter `gammaA`
#' (power function: `gammaA = 1` is linear spending, `gammaA = 2` quadratic)
#' - Hwang, Shi & DeCani $\alpha$-spending (`typeOfDesign = "asHSD"`) with
#' parameter `gammaA` (for details, see Wassmer & Brannath 2016, p. 76)
#' - Standard Wang & Tsiatis Delta classes (`typeOfDesign = "WT"`) and
#' (`typeOfDesign = "WToptim"`)
## ---- include = TRUE, echo = TRUE---------------------------------------------
# Quadratic Kim & DeMets alpha-spending
design <- getDesignGroupSequential(sided = 1, alpha = 0.025,
informationRates = c(0.33, 0.67, 1), typeOfDesign = "asKD", gammaA = 2)
#'
#' ## User-defined $\alpha$-spending functions
#' User-defined $\alpha$-spending functions (`typeOfDesign = "asUser"`) can be
#' obtained via the argument `userAlphaSpending` which must contain a numeric vector
#' with elements $0< \alpha_1 < \ldots < \alpha_{kMax} = \alpha$ that define the
#' values of the cumulative alpha-spending function at each interim analysis.
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
# Example: User-defined alpha-spending function which is very conservative at
# first interim (spend alpha = 0.001), conservative at second (spend an additional
# alpha = 0.01, i.e., total cumulative alpha spent is 0.011 up to second interim),
# and spends the remaining alpha at the final analysis (i.e., cumulative
# alpha = 0.025)
design <- getDesignGroupSequential(sided = 1, alpha = 0.025,
informationRates = c(0.33, 0.67, 1),
typeOfDesign = "asUser",
userAlphaSpending = c(0.001, 0.01 + 0.001, 0.025))
# $stageLevels below extract local significance levels across interim analyses.
# Note that the local significance level is exactly 0.001 at the first
# interim, but slightly >0.01 at the second interim because the design
# exploits correlations between interim analyses.
design$stageLevels
#'
#' # Designs with efficacy and futility interim analyses
#'
#' ## Futility boundaries at interims manually defined on $z$-scale
#' - The argument `futilityBounds` contains a vector of futility bounds (on the
#' $z$-value scale) for each interim (but not the final analysis).
#' - A futility bound of $z = 0$ corresponds to an estimated treatment effect of
#' zero or "null", i.e., in this case futility stopping is recommended if the
#' treatment effect estimate at the interim analysis is zero or "goes in the
#' wrong direction". Futility bounds of $z = -\infty$ (which are numerically
#' equivalent to $z = -6$) correspond to no futility stopping at an interim.
#' - Due to the design of rpact, it is not possible to directly define futility
#' boundaries on the treatment effect scale. If this is desired, one would need
#' to manually convert the treatment effect scale to the $z$-scale or,
#' alternatively, experiment by varying the boundaries on the $z$-scale until
#' this implies the targeted critical values on the treatment effect scale.
#' (Critical values on treatment effect scales are routinely provided by sample
#' size functions for different endpoint types such as `getSampleSizeMeans()` (for
#' continuous endpoints), `getSampleSizeRates()` (for binary endpoints), and
#' `getSampleSizeSurvival()` (for survival endpoints). Please see the R Markdown
#' files for these endpoint types for further details.)
#' - By default, all futility boundaries are non-binding (`bindingFutility = FALSE`).
#' Binding futility boundaries (`bindingFutility = TRUE`) are not recommended.
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
# Example: non-binding futility boundary at each interim in case
# estimated treatment effect is null or goes in "the wrong direction"
design <- getDesignGroupSequential(sided = 1, alpha = 0.025,
informationRates = c(0.33, 0.67, 1), typeOfDesign = "asOF",
futilityBounds = c(0,0), bindingFutility = FALSE)
#'
#' ## Formal $\beta$-spending
#' Formal $\beta$-spending functions are defined in the same way as for
#' $\alpha$-spending functions, e.g., a Pocock type $\beta$-spending can be specified
#' as `typeBetaSpending = "bsP"` and `beta` needs to be specified, the default is
#' `beta = 0.20`. Note that in the case, due to computational constraints, the
#' futility bounds are always non-binding (`bindingFutility = FALSE` produces an
#' error message).
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
# Example: non-binding futility boundary at each interim in case
# estimated treatment effect is null or goes in "the wrong direction"
design <- getDesignGroupSequential(sided = 1, alpha = 0.025, beta = 0.1,
informationRates = c(0.33, 0.67, 1), typeOfDesign = "asOF",
typeBetaSpending = "bsP")
#'
#'
#' # Designs with futility interim analyses only
#' Such designs can be implemented by using a user-defined $\alpha$-spending
#' function which spends all of the Type I error at the final analysis. Note that
#' such designs do not allow stopping for efficacy regardless how persuasive the
#' effect is.
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
# Example: non-binding futility boundary using an O'Brien & Fleming type
# beta spending function. No early stopping for efficacy (i.e., all alpha
# is spent at the final analysis).
design <- getDesignGroupSequential(sided = 1, alpha = 0.025, beta = 0.2,
informationRates = c(0.33, 0.67, 1), typeOfDesign = "asUser",
userAlphaSpending = c(0, 0, 0.025), typeBetaSpending = "bsOF",
bindingFutility = FALSE)
#'
#' # Interpreting, accessing and printing rpact boundary objects
#'
#' ## Information contained in the `design` object
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
design <- getDesignGroupSequential(sided = 1, alpha = 0.025, beta = 0.2,
informationRates = c(0.33, 0.67, 1), typeOfDesign = "asOF",
futilityBounds = c(0, 0), bindingFutility = FALSE)
design
#' The key information is contained in the object including **critical values on the
#' $z$-scale** ("Critical values" in rpact output,
#' `design$criticalValues`) and **one-sided local significance levels** ("Stage levels" in rpact output,`design$stageLevels`).
#' Note that the local significance levels are always given as one-sided levels in
#' rpact even if a two-sided design is specified.
#'
#' `names(design)` provides names of all objects included in the `design` object
#' and `as.data.frame(design)` collects all design information into one data frame.
#' `summary(design)` gives a slightly more detailed output. For more details about
#' applying R generics to rpact objects, please refer to the separte R Markdown
#' file dedicated to this topic.
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
names(design)
#' `summary()` creates a nice presentation of the design:
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
summary(design)
#'
#'
#' ## Stopping probabilities and expected sample size reduction
#' `getDesignCharacteristics(design)` provides more detailed information about the
#' design:
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
designChar <- getDesignCharacteristics(design)
designChar
names(designChar)
#'
#' **Note that the design characteristics depend on beta that needs to be specified
#' in `getDesignGroupSequential`. By default, `beta = 0.20`.**
#'
#' Explanations regarding the output:
#'
#' - **Maximum sample size inflation factor** (`$inflationFactor`): This is the
#' maximal sample size a group-sequential trial requires relative to the sample
#' size of a fixed design without interim analyses.
#' - Probabilities of stopping due to a significant result at each interim or the final analysis
#' (`$rejectionProbabilities`), cumulative power (`$power`), and probability of
#' stopping for futility at each interim (`$futilityProbabilities`). All of these are calculated under the alternative H1.
#' - **Expected sample size** of group-sequential design (relative to fixed design) under
#' the alternative hypothesis H1 (`$averageSampleNumber1`), under the null hypothesis H0
#' (`$averageSampleNumber0`), and under the parameter in the middle between H0 and H1.
#' - In addition, `getDesignCharacteristics(design)` provides the required sample size
#' for an abstract group-sequential single arm trial with a normal outcome, effect size 1,
#' and standard deviation 1 (i.e., the simplest group-sequential setting from a
#' mathematical point of view). The sample size for such a trial without interim
#' analyses is given as
#' `$nFixed` and the maximum sample size of the corresponding group-sequential design as `$shift`.
#'
#'
#' The practical relevance of this abstract design is that the **properties of the
#' design** (critical values, sample size inflation factor, rejection probabilies,
#' etc) **carry over to group-sequential designs regardless of the endpoint (e.g.
#' continuous, binary, or survial)** as they all share the same underlying
#' canonical multivariate normal distribution of the $z$-scores.
#'
#' **Overall stopping probabilities, rejection probabilities, and futility probabilities
#' under the null (H0) and the alternative (H1)** (overall and at each stage) can
#' be calculated using the function `getPowerAndAverageSampleNumber()`. To get
#' these numbers, one needs to provide the maximum sample size and the effect size
#' (0 under H0, 1 under H1) of the corresponding type of design.
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
# theta = 0 for calculations under H0
getPowerAndAverageSampleNumber(design,theta = c(0),nMax = designChar$shift)
# theta = 1 for calculations under alternative H1
getPowerAndAverageSampleNumber(design,theta = 1,nMax = designChar$shift)
#' Note that the power under H0, i.e., the significance level, is slightly below
#' 0.025 in this example as it is calculated under the assumption that the
#' non-binding futility boundaries are adhered to.
#'
#' Both (and even more) values can be obtained with one command
#' `getPowerAndAverageSampleNumber(design,theta = c(0,1),nMax = designChar$shift)`
#'
#' # Plotting rpact boundary objects
#' Boundaries can be plotted using the `plot` (or `plot.TrialDesign`) function
#' which produces a [ggplot2](https://ggplot2.tidyverse.org) object.
#'
#' The most relevant plots for (abstract) boundaries without easily interpretable
#' treatment effect are boundary plots on $z$-scale (`type = 1`) or $p$-value-scale
#' (`type = 3`) as well as plots of the $\alpha$-spending function (`type = 4`).
#' Conveniently, argument `showSource = TRUE`also provides the source data for the
#' plot. For examples of all available plots, see the R Markdown document [How to
#' create admirable plots with
#' rpact](https://vignettes.rpact.org/html/rpact_plot_examples.html).
#'
## ---- include = TRUE, echo = TRUE---------------------------------------------
plot(design, type = 1, showSource = TRUE)
plot(design, type = 3, showSource = TRUE)
plot(design, type = 4, showSource = TRUE)
#'
#' # Comparison of multiple designs
#' Multiple designs can be combined into a design set (`getDesignSet()`) and their
#' properties plotted jointly:
## ---- include = TRUE, echo = TRUE---------------------------------------------
# O'Brien & Fleming, 3 equally spaced stages
d1 <- getDesignGroupSequential(typeOfDesign = "OF", kMax = 3)
# Pocock
d2 <- getDesignGroupSequential(typeOfDesign = "P", kMax = 3)
designSet <- getDesignSet(designs = c(d1, d2), variedParameters = "typeOfDesign")
plot(designSet, type = 1)
#' Even simpler, in rpact 3.0, you can also use `plot(d1, d2)`.
#'
#' ***
#'
#' System: rpact `r packageVersion("rpact")`, `r R.version.string`, platform: `r R.version$platform`
#'
## ---- include = TRUE, echo = FALSE, results = 'asis'--------------------------
printCitation()
#'
#'