Nothing
R/f_design_group_sequential.R
In rpact: Confirmatory Adaptive Clinical Trial Design and Analysis
Defines functions
getSimulatedRejectionsDelayedResponse
.getSimulatedRejectionsDelayedResponse
getPowerAndAverageSampleNumber
.getAverageSampleNumber
.getDesignCharacteristics
getDesignCharacteristics
.getFixedSampleSize
getDesignGroupSequential
.assertIsValidBetaSpent
.assertIsValidAlphaSpent
.getDesignGroupSequential
.getDesignInverseNormalDefaultValues
.getDesignGroupSequentialDefaultValues
.getDesignInverseNormal
getDesignInverseNormal
.calculateDecisionProbabilities
.getDesignGroupSequentialUserDefinedBetaSpending
.getDesignGroupSequentialBetaSpending
.fillVec
.getDesignGroupSequentialBetaSpendingApproaches
.getDesignGroupSequentialUserDefinedAlphaSpending
.getDesignGroupSequentialAlphaSpending
.getDesignGroupSequentialWangAndTsiatisOptimum
.getOptimumDesign
.getDesignGroupSequentialHaybittleAndPeto
.calculateAlphaSpent
.getDesignGroupSequentialPampallonaTsiatis
.getDesignGroupSequentialWangAndTsiatis
.getDesignGroupSequentialKMax1
.createDesign
.validateBaseParameters
.validateTypeOfDesign
.validateGroupSequentialProbabilityResults
.validateGroupSequentialProbabilityResultsMulti
getGroupSequentialProbabilities
.getGroupSequentialProbabilities
Documented in
getDesignCharacteristics
getDesignGroupSequential
getDesignInverseNormal
getGroupSequentialProbabilities
getPowerAndAverageSampleNumber
## |
## | *Group sequential design*
## |
## | This file is part of the R package rpact:
## | Confirmatory Adaptive Clinical Trial Design and Analysis
## |
## | Author: Gernot Wassmer, PhD, and Friedrich Pahlke, PhD
## | Licensed under "GNU Lesser General Public License" version 3
## | License text can be found here: https://www.r-project.org/Licenses/LGPL-3
## |
## | RPACT company website: https://www.rpact.com
## | rpact package website: https://www.rpact.org
## |
## | Contact us for information about our services: info@rpact.com
## |
## | File version: $Revision: 7139 $
## | Last changed: $Date: 2023-06-28 08:15:31 +0200 (Mi, 28 Jun 2023) $
## | Last changed by: $Author: pahlke $
## |
#' @include f_core_constants.R
NULL
.getGroupSequentialProbabilities <- function(decisionMatrix, informationRates) {
return(getGroupSequentialProbabilitiesCpp(decisionMatrix, informationRates))
}
#' @title
#' Get Group Sequential Probabilities
#'
#' @description
#' Calculates probabilities in the group sequential setting.
#'
#' @param decisionMatrix A matrix with either 2 or 4 rows and kMax = length(informationRates) columns, see details.
#' @inheritParams param_informationRates
#'
#' @details
#' Given a sequence of information rates (fixing the correlation structure), and
#' decisionMatrix with either 2 or 4 rows and kMax = length(informationRates) columns,
#' this function calculates a probability matrix containing, for two rows, the probabilities:\cr
#' P(Z_1 <- l_1), P(l_1 <- Z_1 < u_1, Z_2 < l_1),..., P(l_kMax-1 <- Z_kMax-1 < u_kMax-1, Z_kMax < l_l_kMax)\cr
#' P(Z_1 <- u_1), P(l_1 <- Z_1 < u_1, Z_2 < u_1),..., P(l_kMax-1 <- Z_kMax-1 < u_kMax-1, Z_kMax < u_l_kMax)\cr
#' P(Z_1 <- Inf), P(l_1 <- Z_1 < u_1, Z_2 < Inf),..., P(l_kMax-1 <- Z_kMax-1 < u_kMax-1, Z_kMax < Inf)\cr
#' with continuation matrix\cr
#' l_1,...,l_kMax\cr
#' u_1,...,u_kMax\cr
#' For 4 rows, the continuation region contains of two regions and the probability matrix is
#' obtained analogously (cf., Wassmer and Brannath, 2016).
#'
#' @family design functions
#'
#' @template examples_get_group_sequential_probabilities
#'
#' @return Returns a numeric matrix containing the probabilities described in the details section.
#'
#' @export
#'
getGroupSequentialProbabilities <- function(decisionMatrix, informationRates) {
.assertAreValidInformationRates(informationRates)
.assertIsValidDecisionMatrix(decisionMatrix, length(informationRates))
return(.getGroupSequentialProbabilities(decisionMatrix = decisionMatrix, informationRates = informationRates))
}
.validateGroupSequentialProbabilityResultsMulti <- function(...) {
args <- list(...)
for (variableName in names(args)) {
if (!.validateGroupSequentialProbabilityResults(results = args[[variableName]], variableName)) {
return(invisible())
}
}
}
.validateGroupSequentialProbabilityResults <- function(results, variableName) {
numberOfNAs <- sum(is.na(results))
if (numberOfNAs == 0) {
return(TRUE)
}
warning(sprintf(
paste0(
C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"in .getGroupSequentialProbabilities(): ",
"variable '%s' contains %s NA's (%.1f%%)"
),
variableName, numberOfNAs, 100 * numberOfNAs / length(results)
), call. = FALSE)
return(FALSE)
}
.validateTypeOfDesign <- function(design) {
.assertDesignParameterExists(design, "typeOfDesign", C_DEFAULT_TYPE_OF_DESIGN)
design$.setParameterType("userAlphaSpending", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("userBetaSpending", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("deltaWT", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("deltaPT1", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("deltaPT0", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("optimizationCriterion", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("gammaA", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("gammaB", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("typeBetaSpending", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("constantBoundsHP", C_PARAM_NOT_APPLICABLE)
if (!(design$typeOfDesign %in% .getDesignTypes())) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"type of design (", design$typeOfDesign, ") must be one of the following: ", .printDesignTypes()
)
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT) {
.assertDesignParameterExists(design, "deltaWT", NA_real_)
.assertIsSingleNumber(design$deltaWT, "deltaWT", naAllowed = FALSE)
.assertIsInClosedInterval(design$deltaWT, "deltaWT", lower = -0.5, upper = 1)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_PT) {
.assertDesignParameterExists(design, "deltaPT1", NA_real_)
.assertIsSingleNumber(design$deltaPT1, "deltaPT1", naAllowed = FALSE)
.assertIsInClosedInterval(design$deltaPT1, "deltaPT1", lower = -0.5, upper = 1)
.assertDesignParameterExists(design, "deltaPT0", NA_real_)
.assertIsSingleNumber(design$deltaPT0, "deltaPT0", naAllowed = FALSE)
.assertIsInClosedInterval(design$deltaPT0, "deltaPT0", lower = -0.5, upper = 1)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT_OPTIMUM) {
.assertDesignParameterExists(design, "optimizationCriterion", C_OPTIMIZATION_CRITERION_DEFAULT)
if (!.isOptimizationCriterion(design$optimizationCriterion)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"optimization criterion must be one of the following: ", .printOptimizationCriterion()
)
}
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_HP) {
.assertDesignParameterExists(design, "constantBoundsHP", C_CONST_BOUND_HP_DEFAULT)
.assertIsSingleNumber(design$constantBoundsHP, "constantBoundsHP")
.assertIsInClosedInterval(design$constantBoundsHP, "constantBoundsHP", lower = 2, upper = NULL)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_AS_KD) {
.assertDesignParameterExists(design, "gammaA", NA_real_)
.assertIsSingleNumber(design$gammaA, "gammaA", naAllowed = FALSE)
if (design$gammaA < 0.4 || design$gammaA > 8) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"parameter 'gammaA' (", design$gammaA, ") for Kim & DeMets alpha ",
"spending function is out of bounds [0.4; 8]"
)
}
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_AS_HSD) {
.assertDesignParameterExists(design, "gammaA", NA_real_)
.assertIsSingleNumber(design$gammaA, "gammaA", naAllowed = FALSE)
if (design$gammaA < -10 || design$gammaA > 5) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"Parameter 'gammaA' (", design$gammaA, ") for Hwang, Shih & DeCani ",
"alpha spending function is out of bounds [-10; 5]"
)
}
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_AS_USER) {
.validateUserAlphaSpending(design)
design$.setParameterType("userAlphaSpending", C_PARAM_USER_DEFINED)
}
if (.isUndefinedArgument(design$alpha)) {
design$alpha <- C_ALPHA_DEFAULT
design$.setParameterType("alpha", C_PARAM_DEFAULT_VALUE)
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY) {
.assertIsValidAlpha(design$alpha)
design$.setParameterType("userAlphaSpending", C_PARAM_DEFAULT_VALUE)
}
if ((.isBetaSpendingDesignType(design$typeBetaSpending) || !.isAlphaSpendingDesignType(design$typeOfDesign)) &&
(design$informationRates[length(design$informationRates)] != 1)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"For specified design, last information rate should be equal 1"
)
}
if (.isAlphaSpendingDesignType(design$typeOfDesign)) {
.assertDesignParameterExists(design, "typeBetaSpending", C_TYPE_OF_DESIGN_BS_NONE)
if (!.isBetaSpendingDesignType(design$typeBetaSpending, noneIncluded = TRUE)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"type of beta spending must be one of the following: ", .printBetaSpendingDesignTypes()
)
}
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_KD) {
.assertDesignParameterExists(design, "gammaB", NA_real_)
.assertIsSingleNumber(design$gammaB, "gammaB", naAllowed = FALSE)
if (design$gammaB < 0.4 || design$gammaB > 8) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"parameter 'gammaB' (", design$gammaB, ") for Kim & DeMets beta ",
"spending function out of bounds [0.4; 8]"
)
}
}
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_HSD) {
.assertDesignParameterExists(design, "gammaB", NA_real_)
.assertIsSingleNumber(design$gammaB, "gammaB", naAllowed = FALSE)
if (design$gammaB < -10 || design$gammaB > 5) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"parameter 'gammaB' (", design$gammaB, ") for Hwang, Shih & DeCani ",
"beta spending out of bounds [-10; 5]"
)
}
}
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_USER) {
.validateUserBetaSpending(design)
design$.setParameterType("userBetaSpending", C_PARAM_USER_DEFINED)
}
} else {
if (.designParameterExists(design, "typeBetaSpending") && design$typeBetaSpending != C_TYPE_OF_DESIGN_BS_NONE) {
warning("'typeBetaSpending' (", design$typeBetaSpending, ") will be ignored ",
"because 'typeOfDesign' (", design$typeOfDesign, ") is not an alpha spending design",
call. = FALSE
)
design$typeBetaSpending <- C_TYPE_OF_DESIGN_BS_NONE
design$.setParameterType("typeBetaSpending", C_PARAM_DEFAULT_VALUE)
}
if (.designParameterExists(design, "userBetaSpending")) {
userBetaSpending <- NA_real_
warning("'userBetaSpending' (", .arrayToString(design$userBetaSpending), ") will be ignored ",
"because 'typeOfDesign' (", design$typeOfDesign, ") is not an alpha spending design",
call. = FALSE
)
}
}
if (.isUndefinedArgument(design$beta)) {
design$beta <- C_BETA_DEFAULT
design$.setParameterType("beta", C_PARAM_DEFAULT_VALUE)
}
return(invisible(design))
}
.validateBaseParameters <- function(design, twoSidedWarningForDefaultValues = TRUE) {
if (.isDefinedArgument(design$kMax)) {
.assertDesignParameterExists(design, "kMax", C_KMAX_DEFAULT)
.assertIsValidKMax(design$kMax)
if (.isDefinedArgument(design$informationRates)) {
.assertAreValidInformationRates(design$informationRates, design$kMax)
}
if (.isDefinedArgument(design$futilityBounds)) {
.assertAreValidFutilityBounds(design$futilityBounds, design$kMax)
}
}
.assertDesignParameterExists(design, "sided", C_SIDED_DEFAULT)
.assertIsValidSidedParameter(design$sided)
.setKmaxBasedOnAlphaSpendingDefintion(design)
design$informationRates <- .getValidatedInformationRates(design)
design$futilityBounds <- .getValidatedFutilityBounds(design,
twoSidedWarningForDefaultValues = twoSidedWarningForDefaultValues
)
.assertDesignParameterExists(design, "tolerance", C_DESIGN_TOLERANCE_DEFAULT)
if (design$tolerance < 1e-10 || design$tolerance > 1e-03) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'tolerance' (", design$tolerance, ") out of bounds [1e-10; 1e-03]"
)
}
return(invisible(design))
}
.createDesign <- function(...,
designClass,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = C_BINDING_FUTILITY_DEFAULT,
constantBoundsHP = C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
betaAdjustment = NA,
delayedInformation = NA_real_,
tolerance = C_DESIGN_TOLERANCE_DEFAULT) {
.assertIsSingleInteger(kMax, "kMax", naAllowed = TRUE, validateType = FALSE)
.assertIsSingleCharacter(typeOfDesign, "typeOfDesign")
if (typeOfDesign == C_TYPE_OF_DESIGN_AS_USER && !any(is.na(userAlphaSpending))) {
if (!is.na(kMax) && kMax != length(userAlphaSpending)) {
stop(sprintf(
paste0(
C_EXCEPTION_TYPE_CONFLICTING_ARGUMENTS,
"length of 'userAlphaSpending' (%s) must be equal to 'kMax' (%s)"
),
length(userAlphaSpending), kMax
))
}
kMax <- length(userAlphaSpending)
if (kMax > 1 && all(userAlphaSpending[1:(kMax - 1)] == 0)) {
message("Changed type of design to ", sQuote(C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY))
typeOfDesign <- C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY
}
}
.assertIsSingleLogical(bindingFutility, "bindingFutility")
.assertIsNumericVector(delayedInformation, "delayedInformation", naAllowed = TRUE)
.assertIsInClosedInterval(delayedInformation, "delayedInformation", lower = 0, upper = NULL, naAllowed = TRUE)
if (designClass == C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL) {
design <- TrialDesignInverseNormal(
kMax = kMax, bindingFutility = bindingFutility,
delayedInformation = delayedInformation
)
} else if (designClass == C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL) {
design <- TrialDesignGroupSequential(
kMax = kMax, bindingFutility = bindingFutility,
delayedInformation = delayedInformation
)
} else {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'designClass' ('", designClass, "') must be '", C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL, "' or ",
"'", C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL, "'"
)
}
.assertIsSingleInteger(sided, "sided", naAllowed = FALSE, validateType = FALSE)
if (!is.integer(sided) && sided %in% c(1, 2)) {
sided <- as.integer(sided)
}
.assertIsSingleCharacter(optimizationCriterion, "optimizationCriterion")
.assertIsSingleCharacter(typeBetaSpending, "typeBetaSpending")
.assertIsSingleLogical(twoSidedPower, "twoSidedPower", naAllowed = TRUE)
.assertIsSingleLogical(betaAdjustment, "betaAdjustment", naAllowed = TRUE)
.assertIsSingleNumber(alpha, "alpha", naAllowed = TRUE)
.assertIsSingleNumber(beta, "beta", naAllowed = TRUE)
.assertIsSingleNumber(deltaWT, "deltaWT", naAllowed = TRUE)
.assertIsSingleNumber(deltaPT1, "deltaPT1", naAllowed = TRUE)
.assertIsSingleNumber(deltaPT0, "deltaPT0", naAllowed = TRUE)
.assertIsSingleNumber(gammaA, "gammaA", naAllowed = TRUE)
.assertIsSingleNumber(gammaB, "gammaB", naAllowed = TRUE)
.assertIsNumericVector(futilityBounds, "futilityBounds", naAllowed = TRUE)
.assertIsNumericVector(informationRates, "informationRates", naAllowed = TRUE)
.assertIsNumericVector(userAlphaSpending, "userAlphaSpending", naAllowed = TRUE)
.assertIsNumericVector(userBetaSpending, "userBetaSpending", naAllowed = TRUE)
design$alpha <- alpha
design$beta <- beta
design$sided <- sided
design$typeOfDesign <- typeOfDesign
design$deltaWT <- deltaWT
design$deltaPT1 <- deltaPT1
design$deltaPT0 <- deltaPT0
design$gammaA <- gammaA
design$gammaB <- gammaB
design$optimizationCriterion <- optimizationCriterion
design$typeBetaSpending <- typeBetaSpending
design$futilityBounds <- futilityBounds
design$informationRates <- informationRates
design$userAlphaSpending <- userAlphaSpending
design$userBetaSpending <- userBetaSpending
design$bindingFutility <- bindingFutility
design$delayedInformation <- delayedInformation
if (!all(is.na(delayedInformation)) && any(delayedInformation > 0)) {
design$.setParameterType("delayedInformation", C_PARAM_USER_DEFINED)
}
if (design$typeOfDesign != C_TYPE_OF_DESIGN_WT_OPTIMUM && optimizationCriterion != C_OPTIMIZATION_CRITERION_DEFAULT) {
warning(
"'optimizationCriterion' (", optimizationCriterion,
") will be ignored because it is only applicable for 'typeOfDesign' = \"", C_TYPE_OF_DESIGN_WT_OPTIMUM, "\"",
call. = FALSE
)
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_HP) {
.assertIsSingleNumber(constantBoundsHP, "constantBoundsHP")
.assertIsInClosedInterval(constantBoundsHP, "constantBoundsHP", lower = 2, upper = NULL)
design$constantBoundsHP <- constantBoundsHP
} else if (constantBoundsHP != C_CONST_BOUND_HP_DEFAULT) {
warning(
"'constantBoundsHP' (", constantBoundsHP,
") will be ignored because it is only applicable for 'typeOfDesign' = \"", C_TYPE_OF_DESIGN_HP, "\"",
call. = FALSE
)
}
if (is.na(twoSidedPower)) {
design$twoSidedPower <- C_TWO_SIDED_POWER_DEFAULT
design$.setParameterType("twoSidedPower", C_PARAM_DEFAULT_VALUE)
} else {
design$twoSidedPower <- twoSidedPower
design$.setParameterType("twoSidedPower", ifelse(
twoSidedPower == C_TWO_SIDED_POWER_DEFAULT, C_PARAM_DEFAULT_VALUE, C_PARAM_USER_DEFINED
))
}
if (design$sided == 2 && grepl("^bs", design$typeBetaSpending)) {
if (is.na(betaAdjustment)) {
design$betaAdjustment <- TRUE
design$.setParameterType("betaAdjustment", C_PARAM_DEFAULT_VALUE)
} else {
design$betaAdjustment <- betaAdjustment
design$.setParameterType("betaAdjustment", ifelse(betaAdjustment, C_PARAM_DEFAULT_VALUE, C_PARAM_USER_DEFINED))
}
} else if (!is.na(betaAdjustment)) {
warning(
"'betaAdjustment' (", betaAdjustment,
") will be ignored because it is only applicable for two-sided beta-spending designs",
call. = FALSE
)
}
design$tolerance <- tolerance
return(design)
}
.getDesignGroupSequentialKMax1 <- function(design) {
design$criticalValues <- .getOneMinusQNorm(design$alpha / design$sided)
design$alphaSpent[1] <- design$alpha
return(invisible(design))
}
#
# Wang and Tsiatis design
#
.getDesignGroupSequentialWangAndTsiatis <- function(design) {
if (design$typeOfDesign == C_TYPE_OF_DESIGN_P) {
design$criticalValues <- getDesignGroupSequentialPocockCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance
)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_OF) {
design$criticalValues <- getDesignGroupSequentialOBrienAndFlemingCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance
)
} else {
design$criticalValues <- getDesignGroupSequentialDeltaWTCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance,
design$deltaWT
)
}
.calculateAlphaSpent(design)
return(invisible(design))
}
.getDesignGroupSequentialPampallonaTsiatis <- function(design) {
cppResult <- getDesignGroupSequentialPampallonaTsiatisCpp(
design$tolerance, design$beta, design$alpha, design$kMax,
design$deltaPT0, design$deltaPT1, design$informationRates,
design$sided, design$bindingFutility
)
futilityBounds <- cppResult$futilityBounds
criticalValues <- cppResult$criticalValues
probs <- cppResult$probs
if (design$sided == 1) {
design$betaSpent <- cumsum(probs[1, ])
design$power <- cumsum(probs[3, ] - probs[2, ])
} else {
design$betaSpent <- cumsum(probs[3, ] - probs[2, ])
if (design$twoSidedPower) {
design$power <- cumsum(probs[5, ] - probs[4, ] + probs[1, ])
} else {
design$power <- cumsum(probs[5, ] - probs[4, ])
}
}
design$.setParameterType("betaSpent", C_PARAM_GENERATED)
design$.setParameterType("power", C_PARAM_GENERATED)
design$futilityBounds <- futilityBounds[1:(design$kMax - 1)]
design$criticalValues <- criticalValues
design$.setParameterType("futilityBounds", C_PARAM_GENERATED)
design$.setParameterType("criticalValues", C_PARAM_GENERATED)
.calculateAlphaSpent(design)
design$futilityBounds[design$futilityBounds == 0] <- NA_real_
.assertIsValidBetaSpent(design)
return(invisible(design))
}
.calculateAlphaSpent <- function(design) {
if (design$sided == 2) {
if (design$bindingFutility) {
futilityBounds <- design$futilityBounds
futilityBounds[is.na(futilityBounds)] <- 0
decisionMatrix <- matrix(c(
-design$criticalValues, -futilityBounds, 0,
futilityBounds, 0, design$criticalValues
), nrow = 4, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(-design$criticalValues, design$criticalValues), nrow = 2, byrow = TRUE)
}
} else {
if (design$bindingFutility) {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues
), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
rep(C_FUTILITY_BOUNDS_DEFAULT, design$kMax),
design$criticalValues
), nrow = 2, byrow = TRUE)
}
}
tryCatch(
{
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
if (design$sided == 1) {
design$alphaSpent <- cumsum(probs[3, ] - probs[2, ])
} else if (nrow(decisionMatrix) == 2) {
design$alphaSpent <- cumsum(probs[3, ] - probs[2, ] + probs[1, ])
} else {
design$alphaSpent <- cumsum(probs[5, ] - probs[4, ] + probs[1, ])
}
if (!is.na(design$alphaSpent[design$kMax])) {
design$alphaSpent[design$kMax] <- floor(design$alphaSpent[design$kMax] * 1e8) / 1e8
}
design$.setParameterType("alphaSpent", C_PARAM_GENERATED)
},
error = function(e) {
warning("Failed to calculate 'alphaSpent': ", e, call. = FALSE)
}
)
}
#'
#' Haybittle & Peto design
#'
#' @noRd
#'
.getDesignGroupSequentialHaybittleAndPeto <- function(design) {
scale <- .getOneDimensionalRoot(
function(scale) {
design$criticalValues <- c(rep(design$constantBoundsHP, design$kMax - 1), scale)
if (design$sided == 2) {
decisionMatrix <- matrix(c(-design$criticalValues, design$criticalValues), nrow = 2, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ] + probs[1, ]) - design$alpha)
} else {
if (design$bindingFutility) {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues
), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
rep(C_FUTILITY_BOUNDS_DEFAULT, design$kMax),
design$criticalValues
), nrow = 2, byrow = TRUE)
}
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ]) - design$alpha)
}
},
lower = 0, upper = 8, tolerance = design$tolerance,
callingFunctionInformation = ".getDesignGroupSequentialHaybittleAndPeto"
)
design$criticalValues <- c(rep(design$constantBoundsHP, design$kMax - 1), scale)
.calculateAlphaSpent(design)
return(invisible(design))
}
.getOptimumDesign <- function(deltaWT, design) {
scale <- .getOneDimensionalRoot(
function(scale) {
criticalValues <- scale * design$informationRates^(deltaWT - 0.5)
if (design$sided == 2) {
decisionMatrix <- (matrix(c(-criticalValues, criticalValues), nrow = 2, byrow = TRUE))
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ] + probs[1, ]) - design$alpha)
} else {
if (design$bindingFutility) {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
criticalValues
), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
rep(C_FUTILITY_BOUNDS_DEFAULT, design$kMax),
criticalValues
), nrow = 2, byrow = TRUE)
}
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ]) - design$alpha)
}
},
lower = 0, upper = 5, tolerance = design$tolerance, callingFunctionInformation = ".getOptimumDesign"
)
design$criticalValues <- scale * design$informationRates^(deltaWT - 0.5)
designCharacteristics <- .getDesignCharacteristics(design = design)
y <- NA_integer_
if (design$optimizationCriterion == C_OPTIMIZATION_CRITERION_ASNH1) {
y <- designCharacteristics$averageSampleNumber1
}
if (design$optimizationCriterion == C_OPTIMIZATION_CRITERION_ASNIFH1) {
y <- designCharacteristics$inflationFactor + designCharacteristics$averageSampleNumber1
}
if (design$optimizationCriterion == C_OPTIMIZATION_CRITERION_ASN_SUM) {
y <- designCharacteristics$averageSampleNumber0 +
designCharacteristics$averageSampleNumber01 + designCharacteristics$averageSampleNumber1
}
return(y)
}
#'
#' Optimum design within Wang and Tsiatis class
#'
#' @noRd
#'
.getDesignGroupSequentialWangAndTsiatisOptimum <- function(design) {
.assertDesignParameterExists(design, "optimizationCriterion", C_OPTIMIZATION_CRITERION_DEFAULT)
.assertIsOptimizationCriterion(design$optimizationCriterion)
optimumDesign <- stats::optimize(.getOptimumDesign,
design = design,
interval = c(0, 1), tol = 0.0001
)
design$deltaWT <- round(optimumDesign$minimum, 3)
design$.setParameterType("deltaWT", C_PARAM_GENERATED)
# Recalculation of design characteristics with rounded design$deltaWT
design$criticalValues <- getDesignGroupSequentialDeltaWTCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance,
design$deltaWT
)
designCharacteristics <- .getDesignCharacteristics(design = design)
design$power <- designCharacteristics$power
design$.setParameterType("power", C_PARAM_GENERATED)
.calculateAlphaSpent(design)
return(invisible(design))
}
#'
#' Alpha spending approaches
#'
#' @noRd
#'
.getDesignGroupSequentialAlphaSpending <- function(design, userFunctionCallEnabled) {
design$criticalValues <- getDesignGroupSequentialAlphaSpendingCpp(
design$kMax,
design$alpha,
design$gammaA,
design$typeOfDesign,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance
)
.calculateAlphaSpent(design)
return(.getDesignGroupSequentialBetaSpendingApproaches(design, userFunctionCallEnabled))
}
#'
#' User defined alpha spending approach
#'
#' @noRd
#'
.getDesignGroupSequentialUserDefinedAlphaSpending <- function(design, userFunctionCallEnabled) {
design$criticalValues <- rep(NA_real_, design$kMax)
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY) {
design$userAlphaSpending <- rep(0, design$kMax)
design$userAlphaSpending[design$kMax] <- design$alpha
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY && !design$bindingFutility) {
design$criticalValues[1:(design$kMax - 1)] <- C_QNORM_THRESHOLD
design$criticalValues[design$kMax] <- .getOneMinusQNorm(design$alpha / design$sided)
} else {
design$criticalValues <- getDesignGroupSequentialUserDefinedAlphaSpendingCpp(
design$kMax, design$userAlphaSpending, design$sided,
design$informationRates, design$bindingFutility, design$futilityBounds, design$tolerance
)
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY) {
design$criticalValues[1:(design$kMax - 1)] <- C_QNORM_THRESHOLD
}
}
.calculateAlphaSpent(design)
return(invisible(.getDesignGroupSequentialBetaSpendingApproaches(design, userFunctionCallEnabled)))
}
#'
#' Only for alpha spending approaches
#'
#' @noRd
#'
.getDesignGroupSequentialBetaSpendingApproaches <- function(design, userFunctionCallEnabled) {
# beta spending approaches (additional to alpha spending)!
if (.isBetaSpendingDesignType(design$typeBetaSpending,
userDefinedBetaSpendingIncluded = FALSE, noneIncluded = FALSE
)) {
.getDesignGroupSequentialBetaSpending(design, userFunctionCallEnabled)
}
# User defined beta spending
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_USER) {
.getDesignGroupSequentialUserDefinedBetaSpending(design)
}
return(invisible(design))
}
.fillVec <- function(vec, n) {
return(c(vec, rep(NA_real_, n - length(vec))))
}
#'
#' Beta spending approaches (additional to alpha spending)
#' Find shift with beta spending such that last critical values coincide
#'
#' @noRd
#'
.getDesignGroupSequentialBetaSpending <- function(design, userFunctionCallEnabled) {
cppResult <- getDesignGroupSequentialBetaSpendingCpp(
design$criticalValues,
design$kMax,
design$userAlphaSpending,
design$userBetaSpending,
design$informationRates,
design$bindingFutility,
design$tolerance,
design$typeOfDesign,
design$typeBetaSpending,
design$gammaA,
design$gammaB,
design$alpha,
design$beta,
design$sided,
design$betaAdjustment,
design$twoSidedPower
)
design$futilityBounds <- cppResult$futilityBounds
design$criticalValues <- cppResult$criticalValues
design$betaSpent <- cppResult$betaSpent
design$power <- cppResult$power
if (design$sided == 2) {
.calculateAlphaSpent(design)
}
design$.setParameterType("betaSpent", C_PARAM_GENERATED)
design$.setParameterType("power", C_PARAM_GENERATED)
design$.setParameterType("futilityBounds", C_PARAM_GENERATED)
return(invisible(design))
}
#'
#' User defined beta spending
#'
#' Find shift with beta spending such that last critical values coincide
#'
#' @noRd
#'
.getDesignGroupSequentialUserDefinedBetaSpending <- function(design) {
if (design$typeBetaSpending != C_TYPE_OF_DESIGN_BS_USER) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'typeBetaSpending' ('", design$typeBetaSpending, "') must be '", C_TYPE_OF_DESIGN_BS_USER, "'"
)
}
cppResult <- getDesignGroupSequentialUserDefinedBetaSpendingCpp(
design$criticalValues,
design$kMax,
design$userAlphaSpending,
design$userBetaSpending,
design$sided,
design$informationRates,
design$bindingFutility,
design$tolerance,
design$typeOfDesign,
design$gammaA,
design$alpha,
design$betaAdjustment,
design$twoSidedPower
)
design$futilityBounds <- cppResult$futilityBounds
design$criticalValues <- cppResult$criticalValues
design$betaSpent <- cppResult$betaSpent
design$power <- cppResult$power
if (design$sided == 2) {
.calculateAlphaSpent(design)
}
design$.setParameterType("betaSpent", C_PARAM_GENERATED)
design$.setParameterType("power", C_PARAM_GENERATED)
design$.setParameterType("futilityBounds", C_PARAM_GENERATED)
return(invisible(design))
}
#'
#' Calculate stopping, rejection and futility probabilities for delayed response design
#'
#' @noRd
#'
.calculateDecisionProbabilities <- function(sqrtShift, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues) {
kMax <- length(informationRates)
power <- numeric(kMax)
futilityProbabilities <- numeric(kMax - 1)
stoppingProbabilities <- numeric(kMax - 1)
rejectionProbabilities <- numeric(kMax)
contRegionLower <- c(contRegionLower, decisionCriticalValues[kMax])
for (stage in 1:(kMax)) {
if (!is.na(delayedInformation[stage]) && delayedInformation[stage] > 0) {
# information rate vector in case of recruitment stop at 'stage'
informationRatesUponDelay <- c(
informationRates[1:stage],
informationRates[stage] + delayedInformation[stage]
)
if (stage == 1) {
probs1 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionUpper[stage] - sqrtShift * sqrt(informationRatesUponDelay[1]),
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[2]),
C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
-C_UPPER_BOUNDS_DEFAULT, decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[2]),
contRegionLower[stage] - sqrtShift * sqrt(informationRatesUponDelay[1]), C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
rejectionProbabilities[stage] <- probs1[2, stage + 1] - probs1[1, stage + 1] +
probs2[2, stage + 1] - probs2[1, stage + 1]
power[stage] <- rejectionProbabilities[stage]
futilityProbabilities[stage] <- probs2[2, stage]
stoppingProbabilities[stage] <- probs2[2, stage] + 1 - probs1[1, stage]
} else if (stage < kMax) {
probs1 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift * sqrt(informationRatesUponDelay[1:(stage - 1)]),
contRegionUpper[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage]),
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage + 1]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRatesUponDelay[1:(stage - 1)]),
C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift *
sqrt(informationRatesUponDelay[1:(stage - 1)]), -C_UPPER_BOUNDS_DEFAULT,
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage + 1]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRatesUponDelay[1:(stage - 1)]),
contRegionLower[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage]), C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
rejectionProbabilities[stage] <- probs1[2, stage + 1] - probs1[1, stage + 1] +
probs2[2, stage + 1] - probs2[1, stage + 1]
power[stage] <- sum(rejectionProbabilities[1:stage])
futilityProbabilities[stage] <- probs2[2, stage]
stoppingProbabilities[stage] <- probs2[2, stage] + probs1[2, stage] - probs1[1, stage]
} else {
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRates[stage]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRates
)
rejectionProbabilities[stage] <- probs[2, stage] - probs[1, stage]
power[stage] <- sum(rejectionProbabilities[1:stage])
}
} else {
if (stage == 1) {
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[stage] - sqrtShift * sqrt(informationRates[stage]),
contRegionUpper[stage] - sqrtShift * sqrt(informationRates[stage])
),
nrow = 2, byrow = TRUE
), informationRates[1]
)
} else {
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
contRegionLower[stage] - sqrtShift * sqrt(informationRates[stage]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
contRegionUpper[stage] - sqrtShift * sqrt(informationRates[stage])
),
nrow = 2, byrow = TRUE
), informationRates[1:stage]
)
}
rejectionProbabilities[stage] <- probs[3, stage] - probs[2, stage]
if (stage < kMax) {
futilityProbabilities[stage] <- probs[1, stage]
stoppingProbabilities[stage] <- futilityProbabilities[stage] + rejectionProbabilities[stage]
}
power[stage] <- sum(rejectionProbabilities[1:stage])
}
}
return(list(
probs = probs,
power = power,
rejectionProbabilities = rejectionProbabilities,
futilityProbabilities = futilityProbabilities,
stoppingProbabilities = stoppingProbabilities
))
}
#'
#' @title
#' Get Design Inverse Normal
#'
#' @description
#' Provides adjusted boundaries and defines a group sequential design for its use in
#' the inverse normal combination test.
#'
#' @inheritParams getDesignGroupSequential
#'
#' @template details_group_sequential_design
#'
#' @template return_object_trial_design
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_design_inverse_normal
#'
#' @export
#'
getDesignInverseNormal <- function(...,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = 1L, # C_SIDED_DEFAULT
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = c("OF", "P", "WT", "PT", "HP", "WToptimum", "asP", "asOF", "asKD", "asHSD", "asUser", "noEarlyEfficacy"), # C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = c("ASNH1", "ASNIFH1", "ASNsum"), # C_OPTIMIZATION_CRITERION_DEFAULT
gammaA = NA_real_,
typeBetaSpending = c("none", "bsP", "bsOF", "bsKD", "bsHSD", "bsUser"), # C_TYPE_OF_DESIGN_BS_NONE
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = NA,
betaAdjustment = NA,
constantBoundsHP = 3, # C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
tolerance = 1e-08 # C_DESIGN_TOLERANCE_DEFAULT
) {
.warnInCaseOfUnknownArguments(functionName = "getDesignInverseNormal", ...)
return(.getDesignGroupSequential(
designClass = C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
betaAdjustment = betaAdjustment,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
tolerance = tolerance,
userFunctionCallEnabled = TRUE
))
}
.getDesignInverseNormal <- function(...,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = NA,
betaAdjustment = NA,
constantBoundsHP = C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
tolerance = C_DESIGN_TOLERANCE_DEFAULT) {
.warnInCaseOfUnknownArguments(functionName = "getDesignInverseNormal", ...)
return(.getDesignGroupSequential(
designClass = C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
betaAdjustment = betaAdjustment,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
tolerance = tolerance,
userFunctionCallEnabled = FALSE
))
}
.getDesignGroupSequentialDefaultValues <- function() {
return(list(
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
twoSidedPower = C_TWO_SIDED_POWER_DEFAULT,
tolerance = C_DESIGN_TOLERANCE_DEFAULT
))
}
.getDesignInverseNormalDefaultValues <- function() {
return(.getDesignGroupSequentialDefaultValues())
}
#
# Param: userFunctionCallEnabled if \code{TRUE}, additional parameter validation methods will be called.
#
.getDesignGroupSequential <- function(...,
designClass = C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = C_BINDING_FUTILITY_DEFAULT,
constantBoundsHP = C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
betaAdjustment = NA,
delayedInformation = NA_real_,
tolerance = C_DESIGN_TOLERANCE_DEFAULT,
userFunctionCallEnabled = FALSE) {
typeOfDesign <- .matchArgument(typeOfDesign, C_DEFAULT_TYPE_OF_DESIGN)
optimizationCriterion <- .matchArgument(optimizationCriterion, C_OPTIMIZATION_CRITERION_DEFAULT)
typeBetaSpending <- .matchArgument(typeBetaSpending, C_TYPE_OF_DESIGN_BS_NONE)
if (.isDefinedArgument(kMax, argumentExistsValidationEnabled = userFunctionCallEnabled)) {
.assertIsValidKMax(kMax, showWarnings = TRUE)
if (!is.integer(kMax)) {
kMax <- as.integer(kMax)
}
}
if (is.na(bindingFutility)) {
bindingFutility <- C_BINDING_FUTILITY_DEFAULT
} else if (userFunctionCallEnabled && typeOfDesign != C_TYPE_OF_DESIGN_PT &&
!(typeBetaSpending == "bsP" || typeBetaSpending == "bsOF" || typeBetaSpending == "bsKD" ||
typeBetaSpending == "bsHSD" || typeBetaSpending == "bsUser") &&
((!is.na(kMax) && kMax == 1) || any(is.na(futilityBounds)) ||
(!any(is.na(futilityBounds)) && all(futilityBounds == C_FUTILITY_BOUNDS_DEFAULT)))) {
warning("'bindingFutility' (", bindingFutility, ") will be ignored", call. = FALSE)
}
design <- .createDesign(
designClass = designClass,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
betaAdjustment = betaAdjustment,
delayedInformation = delayedInformation,
tolerance = tolerance
)
if (userFunctionCallEnabled) {
.validateBaseParameters(design, twoSidedWarningForDefaultValues = FALSE)
.validateTypeOfDesign(design)
.assertIsValidTolerance(tolerance)
.assertDesignParameterExists(design, "alpha", C_ALPHA_DEFAULT)
.assertDesignParameterExists(design, "beta", C_BETA_DEFAULT)
.assertDesignParameterExists(design, "sided", C_SIDED_DEFAULT)
.assertDesignParameterExists(design, "typeOfDesign", C_DEFAULT_TYPE_OF_DESIGN)
.assertDesignParameterExists(design, "bindingFutility", C_BINDING_FUTILITY_DEFAULT)
.assertDesignParameterExists(design, "tolerance", C_DESIGN_TOLERANCE_DEFAULT)
if (typeOfDesign != C_TYPE_OF_DESIGN_PT) {
if (!is.na(deltaPT1)) {
warning("'deltaPT1' (", deltaPT1, ") will be ignored", call. = FALSE)
}
if (!is.na(deltaPT0)) {
warning("'deltaPT0' (", deltaPT0, ") will be ignored", call. = FALSE)
}
}
if (typeOfDesign != C_TYPE_OF_DESIGN_WT && !is.na(deltaWT)) {
warning("'deltaWT' (", deltaWT, ") will be ignored", call. = FALSE)
}
if (typeOfDesign != C_TYPE_OF_DESIGN_AS_KD &&
typeOfDesign != C_TYPE_OF_DESIGN_AS_HSD && !is.na(gammaA)) {
warning("'gammaA' (", gammaA, ") will be ignored", call. = FALSE)
}
if (typeBetaSpending != C_TYPE_OF_DESIGN_BS_KD &&
typeBetaSpending != C_TYPE_OF_DESIGN_BS_HSD && !is.na(gammaB)) {
warning("'gammaB' (", gammaB, ") will be ignored", call. = FALSE)
}
if (typeBetaSpending != C_TYPE_OF_DESIGN_BS_USER && !all(is.na(userBetaSpending))) {
warning("'userBetaSpending' (", .arrayToString(userBetaSpending), ") will be ignored", call. = FALSE)
}
if (!(typeOfDesign %in% c(C_TYPE_OF_DESIGN_AS_USER)) &&
!all(is.na(userAlphaSpending))) {
warning("'userAlphaSpending' (", .arrayToString(userAlphaSpending), ") will be ignored", call. = FALSE)
}
}
if (design$sided == 2 && design$bindingFutility && design$typeOfDesign != C_TYPE_OF_DESIGN_PT &&
!.isBetaSpendingDesignType(design$typeBetaSpending)) {
warning("'bindingFutility' will be ignored because the test is defined as two-sided", call. = FALSE)
design$bindingFutility <- FALSE
}
if (design$sided == 1 && design$twoSidedPower) {
warning("'twoSidedPower' will be ignored because the test is defined as one-sided", call. = FALSE)
design$twoSidedPower <- FALSE
}
if (userFunctionCallEnabled) {
.validateAlphaAndBeta(design)
}
design$alphaSpent <- rep(NA_real_, design$kMax)
design$betaSpent <- rep(NA_real_, design$kMax)
design$power <- rep(NA_real_, design$kMax)
if (userFunctionCallEnabled) {
design$.setParameterType("betaSpent", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("power", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("alphaSpent", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("criticalValues", C_PARAM_GENERATED)
}
if (design$kMax == 1) {
.getDesignGroupSequentialKMax1(design)
} else {
# Wang and Tsiatis design
if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT ||
design$typeOfDesign == C_TYPE_OF_DESIGN_P ||
design$typeOfDesign == C_TYPE_OF_DESIGN_OF) {
.getDesignGroupSequentialWangAndTsiatis(design)
}
# Pampallona & Tsiatis design
else if (design$typeOfDesign == C_TYPE_OF_DESIGN_PT) {
.getDesignGroupSequentialPampallonaTsiatis(design)
}
# Haybittle & Peto design
else if (design$typeOfDesign == C_TYPE_OF_DESIGN_HP) {
.getDesignGroupSequentialHaybittleAndPeto(design)
}
# Optimum design within Wang and Tsiatis class
else if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT_OPTIMUM) {
.getDesignGroupSequentialWangAndTsiatisOptimum(design)
}
# alpha spending approaches
else if (.isAlphaSpendingDesignType(design$typeOfDesign, userDefinedAlphaSpendingIncluded = FALSE)) {
.getDesignGroupSequentialAlphaSpending(design, userFunctionCallEnabled)
}
# user defined alpha spending approach
else if (design$typeOfDesign %in% c(C_TYPE_OF_DESIGN_AS_USER, C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY)) {
.getDesignGroupSequentialUserDefinedAlphaSpending(design, userFunctionCallEnabled)
} else {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "no calculation routine defined for ", design$typeOfDesign)
}
}
design$stageLevels <- 1 - stats::pnorm(design$criticalValues)
design$.setParameterType("stageLevels", C_PARAM_GENERATED)
if (design$kMax == 1) {
design$.setParameterType("futilityBounds", C_PARAM_NOT_APPLICABLE)
}
if (!all(is.na(design$futilityBounds))) {
if (length(design$futilityBounds) == 0 || all(design$futilityBounds == C_FUTILITY_BOUNDS_DEFAULT)) {
design$.setParameterType("bindingFutility", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("futilityBounds", C_PARAM_NOT_APPLICABLE)
} else if (userFunctionCallEnabled &&
any(design$futilityBounds > design$criticalValues[1:(design$kMax - 1)] - 0.01, na.rm = TRUE)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'futilityBounds' (", .arrayToString(design$futilityBounds), ") too extreme for this situation"
)
}
}
.assertIsValidAlphaSpent(design, userFunctionCallEnabled)
design$.initStages()
# we use 7.5 instead of C_QNORM_THRESHOLD as threshold
design$criticalValues[!is.na(design$criticalValues) & design$criticalValues <= -7.5] <- -Inf
design$criticalValues[!is.na(design$criticalValues) & design$criticalValues >= 7.5] <- Inf
design$futilityBounds[!is.na(design$futilityBounds) & design$futilityBounds <=
C_FUTILITY_BOUNDS_DEFAULT] <- C_FUTILITY_BOUNDS_DEFAULT
if (design$kMax == 1) {
if (!identical(design$informationRates, 1)) {
if (!is.na(design$informationRates)) {
warning("Information rate", ifelse(length(design$informationRates) != 1, "s", ""), " ",
.arrayToString(design$informationRates, vectorLookAndFeelEnabled = TRUE),
" will be ignored",
call. = FALSE
)
}
design$informationRates <- 1
}
design$.setParameterType("informationRates", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("stages", C_PARAM_NOT_APPLICABLE)
}
.assertIsNumericVector(delayedInformation, "delayedInformation", naAllowed = TRUE)
if (all(is.na(delayedInformation))) {
# delayed response design is disabled
return(design)
}
if (all(!is.na(delayedInformation)) && all(delayedInformation < 1e-03)) {
warning("At least one delayed information value must be >= 1e-03 to enable delayed response.",
" 'delayedInformation' (", .arrayToString(delayedInformation), ") will be ignored",
call. = FALSE
)
return(design)
}
# proceed with delayed response design
.assertIsInClosedInterval(delayedInformation, "delayedInformation", lower = 0, upper = NULL)
kMax <- design$kMax
contRegionUpper <- design$criticalValues
contRegionLower <- design$futilityBounds
informationRates <- design$informationRates
decisionCriticalValues <- numeric(kMax)
reversalProbabilities <- numeric(kMax - 1)
if (!all(is.na(delayedInformation)) &&
(design$sided != 1 || all(design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT + 1e-06))) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"decision critical values for delayed response design are only available for ",
"one-sided designs with valid futility bounds"
)
}
if (length(delayedInformation) == 1) {
delayedInformation <- rep(delayedInformation, kMax - 1)
}
if (length(delayedInformation) != kMax - 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'delayedInformation' (", .arrayToString(delayedInformation), ") must have length ",
(kMax - 1), " (kMax - 1)"
)
}
indices <- which(delayedInformation > 0 & delayedInformation < 1e-03)
n <- length(indices)
if (n > 0) {
warning("The", ifelse(n == 1, "", paste0(" ", n)), " delayed information value", ifelse(n == 1, "", "s"), " ",
.arrayToString(delayedInformation[indices], mode = "and"),
" will be replaced by 1e-03 to achieve reasonable results",
call. = FALSE
)
delayedInformation[indices] <- 1e-03
}
# sensible interim choices are restricted by amount of delayed information
eps <- design$informationRates[1:(design$kMax - 1)] + delayedInformation
if (!any(is.na(eps)) && any(eps >= 1)) {
stages <- which(eps >= 1)
stagesInfo <- ifelse(length(stages) == 1, paste0(" ", stages), paste0("s ", .arrayToString(stages, mode = "and")))
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'delayedInformation[stage] + informationRates[stage]' for ",
"stage", stagesInfo, " too large (>= 1). Recruitment stop analysis information + pipeline data ",
"information cannot exceed overall trial information. Instead, the ",
"recruitment stop analysis would be skipped, directly proceeding to the ",
"final analysis"
)
}
# loop iterating through the stages calculation the decision critical values
for (stage in 1:(kMax - 1)) {
if (!is.na(delayedInformation[stage]) && delayedInformation[stage] >= 1e-03 - 1e-06) {
# information rate vector in case of recruitment stop at 'stage'
informationRatesUponDelay <- c(
informationRates[1:stage],
informationRates[stage] + delayedInformation[stage]
)
if (stage == 1) {
decisionCriticalValues[stage] <- .getOneDimensionalRoot(function(secondCriticalValue) {
probs1 <- .getGroupSequentialProbabilities(
matrix(c(contRegionUpper[stage], secondCriticalValue, C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
-C_UPPER_BOUNDS_DEFAULT, secondCriticalValue,
contRegionLower[stage], C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
return(probs1[1, stage + 1] - probs2[2, stage + 1] + probs2[1, stage + 1])
}, lower = -C_UPPER_BOUNDS_DEFAULT, upper = C_UPPER_BOUNDS_DEFAULT, tolerance = design$tolerance)
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionUpper[stage], decisionCriticalValues[stage],
C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
} else {
decisionCriticalValues[stage] <- .getOneDimensionalRoot(function(secondCriticalValue) {
probs1 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)], contRegionUpper[stage], secondCriticalValue,
contRegionUpper[1:(stage - 1)], C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)], -C_UPPER_BOUNDS_DEFAULT, secondCriticalValue,
contRegionUpper[1:(stage - 1)], contRegionLower[stage], C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
return(probs1[1, stage + 1] - probs2[2, stage + 1] + probs2[1, stage + 1])
}, lower = -C_UPPER_BOUNDS_DEFAULT, upper = C_UPPER_BOUNDS_DEFAULT, tolerance = design$tolerance)
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)], contRegionUpper[stage], decisionCriticalValues[stage],
contRegionUpper[1:(stage - 1)], C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
}
if (stage < kMax) {
reversalProbabilities[stage] <- probs[1, stage + 1]
}
} else {
decisionCriticalValues[stage] <- NA_real_
reversalProbabilities[stage] <- NA_real_
}
decisionCriticalValues[kMax] <- contRegionUpper[kMax]
alphaSpent <- .calculateDecisionProbabilities(
sqrtShift = 0, informationRates,
delayedInformation, contRegionUpper, contRegionLower, decisionCriticalValues
)$power
}
decisionCriticalValues[decisionCriticalValues <= -C_UPPER_BOUNDS_DEFAULT + 1e-06] <- NA_real_
decisionCriticalValues[decisionCriticalValues >= C_UPPER_BOUNDS_DEFAULT - 1e-06] <- NA_real_
design$decisionCriticalValues <- decisionCriticalValues
design$reversalProbabilities <- reversalProbabilities
design$delayedInformation <- delayedInformation
design$delayedInformation[design$delayedInformation < 1e-03] <- 0
design$.setParameterType("decisionCriticalValues", C_PARAM_GENERATED)
design$.setParameterType("reversalProbabilities", C_PARAM_GENERATED)
warning("The delayed information design feature is experimental and ",
"hence not fully validated (see www.rpact.com/experimental)",
call. = FALSE
)
return(design)
}
# to avoid error messages in case of repeated p-values computation
.assertIsValidAlphaSpent <- function(design, userFunctionCallEnabled = TRUE) {
if (!userFunctionCallEnabled) {
return(invisible())
}
if (design$informationRates[design$kMax] != 1) {
return(invisible())
}
if (is.na(design$alphaSpent[design$kMax]) || abs(design$alphaSpent[design$kMax] - design$alpha) > 1e-05) {
stop(
C_EXCEPTION_TYPE_RUNTIME_ISSUE, "critical values cannot be calculated ",
"(alpha: ", design$alpha, "; alpha spent at maximum stage: ", design$alphaSpent[design$kMax], ")"
)
}
}
.assertIsValidBetaSpent <- function(design, ..., userFunctionCallEnabled = TRUE, iteration = 1) {
if (!userFunctionCallEnabled) {
return(invisible())
}
if (design$informationRates[design$kMax] != 1) {
return(invisible())
}
if (iteration < 0) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "critical values cannot be calculated")
}
if (is.na(design$betaSpent[design$kMax]) || abs(design$betaSpent[design$kMax] - design$beta) > 1e-05) {
stop(
C_EXCEPTION_TYPE_RUNTIME_ISSUE, "critical values cannot be calculated ",
"(beta spent at maximum stage: ", design$betaSpent[design$kMax], ")"
)
}
}
#'
#' @title
#' Get Design Group Sequential
#'
#' @description
#' Provides adjusted boundaries and defines a group sequential design.
#'
#' @inheritParams param_kMax
#' @inheritParams param_alpha
#' @inheritParams param_beta
#' @inheritParams param_sided
#' @inheritParams param_typeOfDesign
#' @inheritParams param_informationRates
#' @param futilityBounds The futility bounds, defined on the test statistic z scale
#' (numeric vector of length \code{kMax - 1}).
#' @inheritParams param_bindingFutility
#' @param deltaWT Delta for Wang & Tsiatis Delta class.
#' @param deltaPT1 Delta1 for Pampallona & Tsiatis class rejecting H0 boundaries.
#' @param deltaPT0 Delta0 for Pampallona & Tsiatis class rejecting H1 boundaries.
#' @param constantBoundsHP The constant bounds up to stage \code{kMax - 1} for the
#' Haybittle & Peto design (default is \code{3}).
#' @param optimizationCriterion Optimization criterion for optimum design within
#' Wang & Tsiatis class (\code{"ASNH1"}, \code{"ASNIFH1"},
#' \code{"ASNsum"}), default is \code{"ASNH1"}, see details.
#' @param typeBetaSpending Type of beta spending. Type of of beta spending is one of the following:
#' O'Brien & Fleming type beta spending, Pocock type beta spending,
#' Kim & DeMets beta spending, Hwang, Shi & DeCani beta spending, user defined
#' beta spending (\code{"bsOF"}, \code{"bsP"}, \code{"bsKD"},
#' \code{"bsHSD"}, \code{"bsUser"}, default is \code{"none"}).
#' @param gammaA Parameter for alpha spending function.
#' @param gammaB Parameter for beta spending function.
#' @inheritParams param_userAlphaSpending
#' @param delayedInformation Delay of information for delayed response designs. Can be a numeric value or a
#' numeric vector of length \code{kMax - 1}
#' @param userBetaSpending The user defined beta spending. Vector of length \code{kMax} containing the cumulative
#' beta-spending up to each interim stage.
#' @param twoSidedPower For two-sided testing, if \code{twoSidedPower = TRUE} is specified
#' the sample size calculation is performed by considering both tails of the distribution.
#' Default is \code{FALSE}, i.e., it is assumed that one tail probability is equal to 0 or the power
#' should be directed to one part.
#' @param betaAdjustment For two-sided beta spending designs, if \code{betaAdjustement = TRUE} a linear
#' adjustment of the beta spending values is performed if an overlapping of decision regions for futility
#' stopping at earlier stages occurs, otherwise no adjustment is performed (default is \code{TRUE}).
#' @param tolerance The numerical tolerance, default is \code{1e-08}.
#' @inheritParams param_three_dots
#'
#' @template details_group_sequential_design
#'
#' @template return_object_trial_design
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_design_group_sequential
#'
#' @export
#'
getDesignGroupSequential <- function(...,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = 1L, # C_SIDED_DEFAULT
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = c("OF", "P", "WT", "PT", "HP", "WToptimum", "asP", "asOF", "asKD", "asHSD", "asUser", "noEarlyEfficacy"), # C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = c("ASNH1", "ASNIFH1", "ASNsum"), # C_OPTIMIZATION_CRITERION_DEFAULT
gammaA = NA_real_,
typeBetaSpending = c("none", "bsP", "bsOF", "bsKD", "bsHSD", "bsUser"), # C_TYPE_OF_DESIGN_BS_NONE
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = NA,
betaAdjustment = NA,
constantBoundsHP = 3, # C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
delayedInformation = NA_real_,
tolerance = 1e-08 # C_DESIGN_TOLERANCE_DEFAULT
) {
.warnInCaseOfUnknownArguments(functionName = "getDesignGroupSequential", ...)
return(.getDesignGroupSequential(
designClass = C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
betaAdjustment = betaAdjustment,
delayedInformation = delayedInformation,
tolerance = tolerance,
userFunctionCallEnabled = TRUE
))
}
.getFixedSampleSize <- function(alpha, beta, sided, twoSidedPower = C_TWO_SIDED_POWER_DEFAULT) {
.assertIsValidAlphaAndBeta(alpha = alpha, beta = beta)
.assertIsValidSidedParameter(sided)
if (sided == 1) {
return((.getOneMinusQNorm(alpha) + .getOneMinusQNorm(beta))^2)
}
if (twoSidedPower) {
n <- .getOneDimensionalRoot(
function(n) {
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n)) -
stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n)) + beta
},
lower = 0,
upper = 2 * (.getOneMinusQNorm(alpha / 2) + .getOneMinusQNorm(beta))^2,
tolerance = 1e-08, callingFunctionInformation = ".getFixedSampleSize"
)
} else {
n <- (.getOneMinusQNorm(alpha / 2) + .getOneMinusQNorm(beta))^2
}
return(n)
}
#' @title
#' Get Design Characteristics
#'
#' @description
#' Calculates the characteristics of a design and returns it.
#'
#' @inheritParams param_design
#' @inheritParams param_three_dots
#'
#' @details
#' Calculates the inflation factor (IF),
#' the expected reduction in sample size under H1, under H0, and under a value in between H0 and H1.
#' Furthermore, absolute information values are calculated
#' under the prototype case testing H0: mu = 0 against H1: mu = 1.
#'
#' @return Returns a \code{\link{TrialDesignCharacteristics}} object.
#' The following generics (R generic functions) are available for this result object:
#' \itemize{
#' \item \code{\link[=names.FieldSet]{names()}} to obtain the field names,
#' \item \code{\link[=print.FieldSet]{print()}} to print the object,
#' \item \code{\link[=summary.ParameterSet]{summary()}} to display a summary of the object,
#' \item \code{\link[=plot.ParameterSet]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.TrialDesignCharacteristics]{as.data.frame()}} to coerce the object to a \code{\link[base]{data.frame}},
#' \item \code{\link[=as.matrix.FieldSet]{as.matrix()}} to coerce the object to a \code{\link[base]{matrix}}.
#' }
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_design_characteristics
#'
#' @export
#'
getDesignCharacteristics <- function(design = NULL, ...) {
if (is.null(design)) {
design <- .getDefaultDesign(..., type = "characteristics")
.warnInCaseOfUnknownArguments(
functionName = "getDesignCharacteristics",
ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = FALSE), ...
)
} else {
.assertIsTrialDesign(design)
.warnInCaseOfUnknownArguments(functionName = "getDesignCharacteristics", ...)
.warnInCaseOfTwoSidedPowerArgument(...)
}
return(.getDesignCharacteristics(design = design, userFunctionCallEnabled = TRUE))
}
.getDesignCharacteristics <- function(..., design, userFunctionCallEnabled = FALSE) {
.assertIsTrialDesignInverseNormalOrGroupSequential(design)
.assertDesignParameterExists(design, "sided", C_SIDED_DEFAULT)
.assertIsValidSidedParameter(design$sided)
if (userFunctionCallEnabled) {
.validateAlphaAndBeta(design = design)
}
design$informationRates <- .getValidatedInformationRates(design, writeToDesign = FALSE)
if ((design$typeOfDesign == C_TYPE_OF_DESIGN_PT ||
.isBetaSpendingDesignType(design$typeBetaSpending)) && design$sided == 2 && design$kMax == 2) {
design$futilityBounds[is.na(design$futilityBounds)] <- 0
}
design$futilityBounds <- .getValidatedFutilityBounds(design,
writeToDesign = FALSE, twoSidedWarningForDefaultValues = FALSE
)
designCharacteristics <- TrialDesignCharacteristics(design = design)
designCharacteristics$rejectionProbabilities <- rep(NA_real_, design$kMax)
designCharacteristics$.setParameterType("rejectionProbabilities", C_PARAM_NOT_APPLICABLE)
designCharacteristics$futilityProbabilities <- rep(NA_real_, design$kMax - 1)
designCharacteristics$.setParameterType("futilityProbabilities", C_PARAM_NOT_APPLICABLE)
nFixed <- .getFixedSampleSize(
alpha = design$alpha, beta = design$beta,
sided = design$sided, twoSidedPower = design$twoSidedPower
)
designCharacteristics$nFixed <- nFixed
designCharacteristics$.setParameterType("nFixed", C_PARAM_GENERATED)
design$criticalValues[design$criticalValues > 7.5] <- 7.5
if (length(design$decisionCriticalValues) > 0) {
design$decisionCriticalValues[!is.na(design$decisionCriticalValues) & design$decisionCriticalValues > 7.5] <- 7.5
}
informationRates <- design$informationRates
if (design$kMax == 1) {
designCharacteristics$shift <- nFixed
designCharacteristics$.setParameterType("shift", C_PARAM_GENERATED)
designCharacteristics$inflationFactor <- designCharacteristics$shift / nFixed
designCharacteristics$.setParameterType("inflationFactor", C_PARAM_GENERATED)
designCharacteristics$power <- 1 - design$beta
designCharacteristics$.setParameterType("power", design$.getParameterType("power"))
designCharacteristics$.setParameterType("information", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType("averageSampleNumber1", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType("averageSampleNumber01", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType("averageSampleNumber0", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType(".probs", C_PARAM_NOT_APPLICABLE)
return(designCharacteristics)
}
if (!any(is.na(design$delayedInformation)) && length(design$decisionCriticalValues) > 0) {
kMax <- design$kMax
contRegionUpper <- design$criticalValues
contRegionLower <- design$futilityBounds
informationRates <- design$informationRates
decisionCriticalValues <- design$decisionCriticalValues
informationRates <- design$informationRates
delayedInformation <- design$delayedInformation
kMax <- length(informationRates)
shift <- .getOneDimensionalRoot(
function(shift) {
resultsH1 <- .calculateDecisionProbabilities(sqrt(shift), informationRates, delayedInformation, contRegionUpper, contRegionLower, decisionCriticalValues)
return(resultsH1$power[kMax] - 1 + design$beta)
},
lower = 0, upper = 4 * nFixed,
tolerance = design$tolerance, callingFunctionInformation = ".getDesignCharacteristics"
)
stopping <- numeric(kMax)
futility <- numeric(kMax)
rejectionProbabilities <- numeric(kMax)
resultsH1 <- .calculateDecisionProbabilities(
sqrtShift = sqrt(shift), informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues
)
resultsH01 <- .calculateDecisionProbabilities(
sqrtShift = sqrt(shift) / 2, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues
)
resultsH0 <- .calculateDecisionProbabilities(
sqrtShift = 0, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues
)
designCharacteristics$shift <- shift
designCharacteristics$.probs <- resultsH1$probs
designCharacteristics$power <- resultsH1$power
designCharacteristics$information <- informationRates * shift
designCharacteristics$averageSampleNumber1 <-
(shift - sum(resultsH1$stoppingProbabilities * (informationRates[kMax] -
delayedInformation - informationRates[1:(kMax - 1)]) * shift)) / nFixed
designCharacteristics$averageSampleNumber01 <-
(shift - sum(resultsH01$stoppingProbabilities * (informationRates[kMax] -
delayedInformation - informationRates[1:(kMax - 1)]) * shift)) / nFixed
designCharacteristics$averageSampleNumber0 <-
(shift - sum(resultsH0$stoppingProbabilities * (informationRates[kMax] -
delayedInformation - informationRates[1:(kMax - 1)]) * shift)) / nFixed
futilityProbabilities <- resultsH1$futilityProbabilities
rejectionProbabilities <- resultsH1$power
stoppingProbabilities <- resultsH1$stoppingProbabilities
rejectionProbabilities[2:kMax] <- resultsH1$power[2:kMax] - rejectionProbabilities[1:(kMax - 1)]
designCharacteristics$rejectionProbabilities <- rejectionProbabilities
designCharacteristics$futilityProbabilities <- futilityProbabilities
} else if ((design$typeOfDesign == C_TYPE_OF_DESIGN_PT ||
.isBetaSpendingDesignType(design$typeBetaSpending)) && design$sided == 2) {
design$futilityBounds[is.na(design$futilityBounds)] <- 0
shift <- .getOneDimensionalRoot(
function(shift) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
c(-design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
c(design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 4, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
if (design$twoSidedPower) {
return(sum(probs[5, ] - probs[4, ] + probs[1, ]) - 1 + design$beta)
} else {
return(sum(probs[5, ] - probs[4, ]) - 1 + design$beta)
}
},
lower = 0, upper = 4 * (.getOneMinusQNorm(design$alpha / design$sided) + .getOneMinusQNorm(design$beta))^2,
tolerance = design$tolerance, callingFunctionInformation = ".getDesignCharacteristics"
)
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
c(-design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
c(design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 4, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$shift <- shift
designCharacteristics$.probs <- probs
if (design$twoSidedPower) {
designCharacteristics$power <- cumsum(probs[5, ] - probs[4, ] + probs[1, ])
designCharacteristics$rejectionProbabilities <- probs[5, ] - probs[4, ] + probs[1, ]
} else {
designCharacteristics$power <- cumsum(probs[5, ] - probs[4, ])
designCharacteristics$rejectionProbabilities <- probs[5, ] - probs[4, ]
}
if (design$kMax > 1) {
designCharacteristics$futilityProbabilities <- probs[3, 1:(design$kMax - 1)] - probs[2, 1:(design$kMax - 1)]
}
designCharacteristics$information <- informationRates * shift
designCharacteristics$averageSampleNumber1 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs, shift, nFixed
)
decisionMatrix <- matrix(c(
-design$criticalValues,
c(-design$futilityBounds, 0),
c(design$futilityBounds, 0),
design$criticalValues
), nrow = 4, byrow = TRUE)
probs0 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber0 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs0, shift, nFixed
)
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates) / 2,
c(-design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]) / 2, 0),
c(design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]) / 2, 0),
design$criticalValues - sqrt(shift * informationRates) / 2
), nrow = 4, byrow = TRUE)
probs01 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber01 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs01, shift, nFixed
)
design$futilityBounds[design$futilityBounds == 0] <- NA_real_
} else {
shift <- .getOneDimensionalRoot(
function(shift) {
if (design$sided == 2) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
if (design$twoSidedPower) {
return(sum(probs[3, ] - probs[2, ] + probs[1, ]) - 1 + design$beta)
} else {
return(sum(probs[3, ] - probs[2, ]) - 1 + design$beta)
}
} else {
shiftedFutilityBounds <- design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)])
shiftedFutilityBounds[design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <-
C_FUTILITY_BOUNDS_DEFAULT
decisionMatrix <- matrix(c(
shiftedFutilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
return(sum(probs[3, ] - probs[2, ]) - 1 + design$beta)
}
},
lower = 0, upper = 4 * (.getOneMinusQNorm(design$alpha / design$sided) + .getOneMinusQNorm(design$beta))^2,
tolerance = design$tolerance, callingFunctionInformation = ".getDesignCharacteristics"
)
if (design$sided == 2) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
} else {
shiftedFutilityBounds <- design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)])
shiftedFutilityBounds[design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <-
C_FUTILITY_BOUNDS_DEFAULT
decisionMatrix <- matrix(c(
shiftedFutilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
}
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$shift <- shift
designCharacteristics$.probs <- probs
if (design$twoSidedPower) {
designCharacteristics$power <- cumsum(probs[3, ] - probs[2, ] + probs[1, ])
designCharacteristics$rejectionProbabilities <- probs[3, ] - probs[2, ] + probs[1, ]
} else {
designCharacteristics$power <- cumsum(probs[3, ] - probs[2, ])
designCharacteristics$rejectionProbabilities <- probs[3, ] - probs[2, ]
}
if (design$kMax > 1) {
if (design$sided == 2) {
designCharacteristics$futilityProbabilities <- rep(0, design$kMax - 1)
} else {
designCharacteristics$futilityProbabilities <- probs[1, 1:(design$kMax - 1)]
designCharacteristics$futilityProbabilities[design$futilityBounds == C_FUTILITY_BOUNDS_DEFAULT] <- 0
}
designCharacteristics$power <- .getNoEarlyEfficacyZeroCorrectedValues(
design, designCharacteristics$power
)
designCharacteristics$rejectionProbabilities <- .getNoEarlyEfficacyZeroCorrectedValues(
design, designCharacteristics$rejectionProbabilities
)
}
designCharacteristics$information <- informationRates * shift
designCharacteristics$averageSampleNumber1 <- .getAverageSampleNumber(
design$kMax,
design$informationRates, probs, shift, nFixed
)
if (design$sided == 2) {
decisionMatrix <- matrix(c(-design$criticalValues, design$criticalValues), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues
), nrow = 2, byrow = TRUE)
}
probs0 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber0 <- .getAverageSampleNumber(
design$kMax,
design$informationRates, probs0, shift, nFixed
)
if (design$sided == 2) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates) / 2,
design$criticalValues - sqrt(shift * informationRates) / 2
), nrow = 2, byrow = TRUE)
} else {
shiftedFutilityBounds <- design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]) / 2
shiftedFutilityBounds[design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <- C_FUTILITY_BOUNDS_DEFAULT
decisionMatrix <- matrix(c(shiftedFutilityBounds, C_FUTILITY_BOUNDS_DEFAULT, design$criticalValues -
sqrt(shift * informationRates) / 2), nrow = 2, byrow = TRUE)
}
probs01 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber01 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs01, shift, nFixed
)
}
design$criticalValues[design$criticalValues >= 7.5 - 1e-8] <- Inf
designCharacteristics$.setParameterType("shift", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("power", C_PARAM_GENERATED)
designCharacteristics$.setParameterType(".probs", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("rejectionProbabilities", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("information", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("futilityProbabilities", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("averageSampleNumber0", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("averageSampleNumber01", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("averageSampleNumber1", C_PARAM_GENERATED)
designCharacteristics$inflationFactor <- shift / nFixed
designCharacteristics$.setParameterType("inflationFactor", C_PARAM_GENERATED)
if (is.na(designCharacteristics$inflationFactor) ||
designCharacteristics$inflationFactor > 4 || designCharacteristics$inflationFactor < 1 - 1e-05) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "inflation factor cannot be calculated")
}
return(designCharacteristics)
}
.getAverageSampleNumber <- function(kMax, informationRates, probs, shift, nFixed) {
if (nrow(probs) == 3) {
return((shift - sum((probs[3, 1:(kMax - 1)] - probs[2, 1:(kMax - 1)] + probs[1, 1:(kMax - 1)]) *
(informationRates[kMax] - informationRates[1:(kMax - 1)]) * shift)) / nFixed)
} else {
return((shift - sum((probs[5, 1:(kMax - 1)] -
probs[4, 1:(kMax - 1)] + probs[3, 1:(kMax - 1)] - probs[2, 1:(kMax - 1)] + probs[1, 1:(kMax - 1)]) *
(informationRates[kMax] - informationRates[1:(kMax - 1)]) * shift)) / nFixed)
}
}
#'
#' @title
#' Get Power And Average Sample Number
#'
#' @description
#' Returns the power and average sample number of the specified design.
#'
#' @inheritParams param_design
#' @inheritParams param_theta
#' @inheritParams param_nMax
#'
#' @details
#' This function returns the power and average sample number (ASN) of the specified
#' design for the prototype case which is testing H0: mu = mu0 in a one-sample design.
#' \code{theta} represents the standardized effect \code{(mu - mu0) / sigma} and power and ASN
#' is calculated for maximum sample size \code{nMax}.
#' For other designs than the one-sample test of a mean the standardized effect needs to be adjusted accordingly.
#'
#' @return Returns a \code{\link{PowerAndAverageSampleNumberResult}} object.
#' The following generics (R generic functions) are available for this result object:
#' \itemize{
#' \item \code{\link[=names.FieldSet]{names()}} to obtain the field names,
#' \item \code{\link[=print.FieldSet]{print()}} to print the object,
#' \item \code{\link[=summary.ParameterSet]{summary()}} to display a summary of the object,
#' \item \code{\link[=plot.ParameterSet]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.PowerAndAverageSampleNumberResult]{as.data.frame()}}
#' to coerce the object to a \code{\link[base]{data.frame}},
#' \item \code{\link[=as.matrix.FieldSet]{as.matrix()}} to coerce the object to a \code{\link[base]{matrix}}.
#' }
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_power_and_average_sample_number
#'
#' @export
#'
getPowerAndAverageSampleNumber <- function(design, theta = seq(-1, 1, 0.02), nMax = 100) {
.assertIsTrialDesign(design)
.assertIsSingleNumber(nMax, "nMax")
.assertIsInClosedInterval(nMax, "nMax", lower = 1, upper = NULL)
return(PowerAndAverageSampleNumberResult(design = design, theta = theta, nMax = nMax))
}
.getSimulatedRejectionsDelayedResponse <- function(delta, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues, iterations, seed = NA_real_) {
seed <- .setSeed(seed)
kMax <- length(informationRates)
zVector <- numeric(kMax)
reject <- 0L
for (i in 1:iterations) {
for (stage in 1:kMax) {
if (stage == 1) {
zVector[stage] <- stats::rnorm(1, delta * sqrt(informationRates[1]), 1)
} else {
zVector[stage] <- (sqrt(informationRates[stage - 1]) * zVector[stage - 1] +
sqrt(informationRates[stage] - informationRates[stage - 1]) *
stats::rnorm(1, delta * sqrt(informationRates[stage] - informationRates[stage - 1]), 1)) /
sqrt(informationRates[stage])
}
if (!is.na(decisionCriticalValues[stage]) && stage < kMax &&
(zVector[stage] > contRegionUpper[stage] || zVector[stage] < contRegionLower[stage])) {
if ((sqrt(informationRates[stage]) * zVector[stage] + sqrt(delayedInformation[stage]) *
stats::rnorm(1, delta * sqrt(delayedInformation[stage]), 1)) /
sqrt(informationRates[stage] + delayedInformation[stage]) > decisionCriticalValues[stage]) {
reject <- reject + 1L
}
break
}
if (stage == kMax && zVector[stage] > decisionCriticalValues[stage]) {
reject <- reject + 1L
}
}
}
simulatedAlpha <- reject / iterations
return(list(
simulatedAlpha = simulatedAlpha,
delta = delta,
iterations = iterations,
seed = seed
))
}
#'
#' @title
#' Simulated Rejections Delayed Response
#'
#' @description
#' Simulates the rejection probability of a delayed response group sequential design with specified parameters.
#'
#' @inheritParams param_design
#' @param delta The delay value.
#' @param iterations The number of simulation iterations.
#' @inheritParams param_seed
#'
#' @details
#' By default, delta = 0, i.e., the Type error rate is simulated.
#'
#' @return Returns a list summarizing the rejection probabilities.
#'
#' @noRd
#'
getSimulatedRejectionsDelayedResponse <- function(design, ..., delta = 0, iterations = 10000, seed = NA_real_) {
.assertIsTrialDesignInverseNormalOrGroupSequential(design)
.assertIsSingleNumber(delta, "delta")
.assertIsValidIterationsAndSeed(iterations, seed, zeroIterationsAllowed = FALSE)
if (!design$.isDelayedResponseDesign()) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'design' must be a delayed response design with specified 'delayedInformation'")
}
startTime <- Sys.time()
result <- .getSimulatedRejectionsDelayedResponse(
delta = delta,
informationRates = design$informationRates,
delayedInformation = design$delayedInformation,
contRegionUpper = design$criticalValues,
contRegionLower = design$futilityBounds,
decisionCriticalValues = design$decisionCriticalValues,
iterations = iterations,
seed = seed
)
simulatedAlpha <- result$simulatedAlpha
stdError <- sqrt(simulatedAlpha * (1 - simulatedAlpha) / iterations)
ciLower <- simulatedAlpha - 1.96 * stdError
ciUpper <- simulatedAlpha + 1.96 * stdError
result$confidenceIntervall <- c(ciLower, ciUpper)
# simulated Type I error rate is within the 95% error bounds
result$alphaWithin95ConfidenceIntervall <- design$alpha > ciLower && design$alpha < ciUpper
result$time <- Sys.time() - startTime
return(result)
}
Try the rpact package in your browser
Any scripts or data that you put into this service are public.
rpact documentation built on July 9, 2023, 6:30 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions getSimulatedRejectionsDelayedResponse .getSimulatedRejectionsDelayedResponse getPowerAndAverageSampleNumber .getAverageSampleNumber .getDesignCharacteristics getDesignCharacteristics .getFixedSampleSize getDesignGroupSequential .assertIsValidBetaSpent .assertIsValidAlphaSpent .getDesignGroupSequential .getDesignInverseNormalDefaultValues .getDesignGroupSequentialDefaultValues .getDesignInverseNormal getDesignInverseNormal .calculateDecisionProbabilities .getDesignGroupSequentialUserDefinedBetaSpending .getDesignGroupSequentialBetaSpending .fillVec .getDesignGroupSequentialBetaSpendingApproaches .getDesignGroupSequentialUserDefinedAlphaSpending .getDesignGroupSequentialAlphaSpending .getDesignGroupSequentialWangAndTsiatisOptimum .getOptimumDesign .getDesignGroupSequentialHaybittleAndPeto .calculateAlphaSpent .getDesignGroupSequentialPampallonaTsiatis .getDesignGroupSequentialWangAndTsiatis .getDesignGroupSequentialKMax1 .createDesign .validateBaseParameters .validateTypeOfDesign .validateGroupSequentialProbabilityResults .validateGroupSequentialProbabilityResultsMulti getGroupSequentialProbabilities .getGroupSequentialProbabilities
Documented in getDesignCharacteristics getDesignGroupSequential getDesignInverseNormal getGroupSequentialProbabilities getPowerAndAverageSampleNumber
## |
## | *Group sequential design*
## |
## | This file is part of the R package rpact:
## | Confirmatory Adaptive Clinical Trial Design and Analysis
## |
## | Author: Gernot Wassmer, PhD, and Friedrich Pahlke, PhD
## | Licensed under "GNU Lesser General Public License" version 3
## | License text can be found here: https://www.r-project.org/Licenses/LGPL-3
## |
## | RPACT company website: https://www.rpact.com
## | rpact package website: https://www.rpact.org
## |
## | Contact us for information about our services: info@rpact.com
## |
## | File version: $Revision: 7139 $
## | Last changed: $Date: 2023-06-28 08:15:31 +0200 (Mi, 28 Jun 2023) $
## | Last changed by: $Author: pahlke $
## |
#' @include f_core_constants.R
NULL
.getGroupSequentialProbabilities <- function(decisionMatrix, informationRates) {
return(getGroupSequentialProbabilitiesCpp(decisionMatrix, informationRates))
}
#' @title
#' Get Group Sequential Probabilities
#'
#' @description
#' Calculates probabilities in the group sequential setting.
#'
#' @param decisionMatrix A matrix with either 2 or 4 rows and kMax = length(informationRates) columns, see details.
#' @inheritParams param_informationRates
#'
#' @details
#' Given a sequence of information rates (fixing the correlation structure), and
#' decisionMatrix with either 2 or 4 rows and kMax = length(informationRates) columns,
#' this function calculates a probability matrix containing, for two rows, the probabilities:\cr
#' P(Z_1 <- l_1), P(l_1 <- Z_1 < u_1, Z_2 < l_1),..., P(l_kMax-1 <- Z_kMax-1 < u_kMax-1, Z_kMax < l_l_kMax)\cr
#' P(Z_1 <- u_1), P(l_1 <- Z_1 < u_1, Z_2 < u_1),..., P(l_kMax-1 <- Z_kMax-1 < u_kMax-1, Z_kMax < u_l_kMax)\cr
#' P(Z_1 <- Inf), P(l_1 <- Z_1 < u_1, Z_2 < Inf),..., P(l_kMax-1 <- Z_kMax-1 < u_kMax-1, Z_kMax < Inf)\cr
#' with continuation matrix\cr
#' l_1,...,l_kMax\cr
#' u_1,...,u_kMax\cr
#' For 4 rows, the continuation region contains of two regions and the probability matrix is
#' obtained analogously (cf., Wassmer and Brannath, 2016).
#'
#' @family design functions
#'
#' @template examples_get_group_sequential_probabilities
#'
#' @return Returns a numeric matrix containing the probabilities described in the details section.
#'
#' @export
#'
getGroupSequentialProbabilities <- function(decisionMatrix, informationRates) {
.assertAreValidInformationRates(informationRates)
.assertIsValidDecisionMatrix(decisionMatrix, length(informationRates))
return(.getGroupSequentialProbabilities(decisionMatrix = decisionMatrix, informationRates = informationRates))
}
.validateGroupSequentialProbabilityResultsMulti <- function(...) {
args <- list(...)
for (variableName in names(args)) {
if (!.validateGroupSequentialProbabilityResults(results = args[[variableName]], variableName)) {
return(invisible())
}
}
}
.validateGroupSequentialProbabilityResults <- function(results, variableName) {
numberOfNAs <- sum(is.na(results))
if (numberOfNAs == 0) {
return(TRUE)
}
warning(sprintf(
paste0(
C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"in .getGroupSequentialProbabilities(): ",
"variable '%s' contains %s NA's (%.1f%%)"
),
variableName, numberOfNAs, 100 * numberOfNAs / length(results)
), call. = FALSE)
return(FALSE)
}
.validateTypeOfDesign <- function(design) {
.assertDesignParameterExists(design, "typeOfDesign", C_DEFAULT_TYPE_OF_DESIGN)
design$.setParameterType("userAlphaSpending", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("userBetaSpending", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("deltaWT", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("deltaPT1", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("deltaPT0", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("optimizationCriterion", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("gammaA", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("gammaB", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("typeBetaSpending", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("constantBoundsHP", C_PARAM_NOT_APPLICABLE)
if (!(design$typeOfDesign %in% .getDesignTypes())) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"type of design (", design$typeOfDesign, ") must be one of the following: ", .printDesignTypes()
)
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT) {
.assertDesignParameterExists(design, "deltaWT", NA_real_)
.assertIsSingleNumber(design$deltaWT, "deltaWT", naAllowed = FALSE)
.assertIsInClosedInterval(design$deltaWT, "deltaWT", lower = -0.5, upper = 1)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_PT) {
.assertDesignParameterExists(design, "deltaPT1", NA_real_)
.assertIsSingleNumber(design$deltaPT1, "deltaPT1", naAllowed = FALSE)
.assertIsInClosedInterval(design$deltaPT1, "deltaPT1", lower = -0.5, upper = 1)
.assertDesignParameterExists(design, "deltaPT0", NA_real_)
.assertIsSingleNumber(design$deltaPT0, "deltaPT0", naAllowed = FALSE)
.assertIsInClosedInterval(design$deltaPT0, "deltaPT0", lower = -0.5, upper = 1)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT_OPTIMUM) {
.assertDesignParameterExists(design, "optimizationCriterion", C_OPTIMIZATION_CRITERION_DEFAULT)
if (!.isOptimizationCriterion(design$optimizationCriterion)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"optimization criterion must be one of the following: ", .printOptimizationCriterion()
)
}
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_HP) {
.assertDesignParameterExists(design, "constantBoundsHP", C_CONST_BOUND_HP_DEFAULT)
.assertIsSingleNumber(design$constantBoundsHP, "constantBoundsHP")
.assertIsInClosedInterval(design$constantBoundsHP, "constantBoundsHP", lower = 2, upper = NULL)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_AS_KD) {
.assertDesignParameterExists(design, "gammaA", NA_real_)
.assertIsSingleNumber(design$gammaA, "gammaA", naAllowed = FALSE)
if (design$gammaA < 0.4 || design$gammaA > 8) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"parameter 'gammaA' (", design$gammaA, ") for Kim & DeMets alpha ",
"spending function is out of bounds [0.4; 8]"
)
}
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_AS_HSD) {
.assertDesignParameterExists(design, "gammaA", NA_real_)
.assertIsSingleNumber(design$gammaA, "gammaA", naAllowed = FALSE)
if (design$gammaA < -10 || design$gammaA > 5) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"Parameter 'gammaA' (", design$gammaA, ") for Hwang, Shih & DeCani ",
"alpha spending function is out of bounds [-10; 5]"
)
}
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_AS_USER) {
.validateUserAlphaSpending(design)
design$.setParameterType("userAlphaSpending", C_PARAM_USER_DEFINED)
}
if (.isUndefinedArgument(design$alpha)) {
design$alpha <- C_ALPHA_DEFAULT
design$.setParameterType("alpha", C_PARAM_DEFAULT_VALUE)
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY) {
.assertIsValidAlpha(design$alpha)
design$.setParameterType("userAlphaSpending", C_PARAM_DEFAULT_VALUE)
}
if ((.isBetaSpendingDesignType(design$typeBetaSpending) || !.isAlphaSpendingDesignType(design$typeOfDesign)) &&
(design$informationRates[length(design$informationRates)] != 1)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"For specified design, last information rate should be equal 1"
)
}
if (.isAlphaSpendingDesignType(design$typeOfDesign)) {
.assertDesignParameterExists(design, "typeBetaSpending", C_TYPE_OF_DESIGN_BS_NONE)
if (!.isBetaSpendingDesignType(design$typeBetaSpending, noneIncluded = TRUE)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"type of beta spending must be one of the following: ", .printBetaSpendingDesignTypes()
)
}
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_KD) {
.assertDesignParameterExists(design, "gammaB", NA_real_)
.assertIsSingleNumber(design$gammaB, "gammaB", naAllowed = FALSE)
if (design$gammaB < 0.4 || design$gammaB > 8) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"parameter 'gammaB' (", design$gammaB, ") for Kim & DeMets beta ",
"spending function out of bounds [0.4; 8]"
)
}
}
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_HSD) {
.assertDesignParameterExists(design, "gammaB", NA_real_)
.assertIsSingleNumber(design$gammaB, "gammaB", naAllowed = FALSE)
if (design$gammaB < -10 || design$gammaB > 5) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"parameter 'gammaB' (", design$gammaB, ") for Hwang, Shih & DeCani ",
"beta spending out of bounds [-10; 5]"
)
}
}
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_USER) {
.validateUserBetaSpending(design)
design$.setParameterType("userBetaSpending", C_PARAM_USER_DEFINED)
}
} else {
if (.designParameterExists(design, "typeBetaSpending") && design$typeBetaSpending != C_TYPE_OF_DESIGN_BS_NONE) {
warning("'typeBetaSpending' (", design$typeBetaSpending, ") will be ignored ",
"because 'typeOfDesign' (", design$typeOfDesign, ") is not an alpha spending design",
call. = FALSE
)
design$typeBetaSpending <- C_TYPE_OF_DESIGN_BS_NONE
design$.setParameterType("typeBetaSpending", C_PARAM_DEFAULT_VALUE)
}
if (.designParameterExists(design, "userBetaSpending")) {
userBetaSpending <- NA_real_
warning("'userBetaSpending' (", .arrayToString(design$userBetaSpending), ") will be ignored ",
"because 'typeOfDesign' (", design$typeOfDesign, ") is not an alpha spending design",
call. = FALSE
)
}
}
if (.isUndefinedArgument(design$beta)) {
design$beta <- C_BETA_DEFAULT
design$.setParameterType("beta", C_PARAM_DEFAULT_VALUE)
}
return(invisible(design))
}
.validateBaseParameters <- function(design, twoSidedWarningForDefaultValues = TRUE) {
if (.isDefinedArgument(design$kMax)) {
.assertDesignParameterExists(design, "kMax", C_KMAX_DEFAULT)
.assertIsValidKMax(design$kMax)
if (.isDefinedArgument(design$informationRates)) {
.assertAreValidInformationRates(design$informationRates, design$kMax)
}
if (.isDefinedArgument(design$futilityBounds)) {
.assertAreValidFutilityBounds(design$futilityBounds, design$kMax)
}
}
.assertDesignParameterExists(design, "sided", C_SIDED_DEFAULT)
.assertIsValidSidedParameter(design$sided)
.setKmaxBasedOnAlphaSpendingDefintion(design)
design$informationRates <- .getValidatedInformationRates(design)
design$futilityBounds <- .getValidatedFutilityBounds(design,
twoSidedWarningForDefaultValues = twoSidedWarningForDefaultValues
)
.assertDesignParameterExists(design, "tolerance", C_DESIGN_TOLERANCE_DEFAULT)
if (design$tolerance < 1e-10 || design$tolerance > 1e-03) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'tolerance' (", design$tolerance, ") out of bounds [1e-10; 1e-03]"
)
}
return(invisible(design))
}
.createDesign <- function(...,
designClass,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = C_BINDING_FUTILITY_DEFAULT,
constantBoundsHP = C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
betaAdjustment = NA,
delayedInformation = NA_real_,
tolerance = C_DESIGN_TOLERANCE_DEFAULT) {
.assertIsSingleInteger(kMax, "kMax", naAllowed = TRUE, validateType = FALSE)
.assertIsSingleCharacter(typeOfDesign, "typeOfDesign")
if (typeOfDesign == C_TYPE_OF_DESIGN_AS_USER && !any(is.na(userAlphaSpending))) {
if (!is.na(kMax) && kMax != length(userAlphaSpending)) {
stop(sprintf(
paste0(
C_EXCEPTION_TYPE_CONFLICTING_ARGUMENTS,
"length of 'userAlphaSpending' (%s) must be equal to 'kMax' (%s)"
),
length(userAlphaSpending), kMax
))
}
kMax <- length(userAlphaSpending)
if (kMax > 1 && all(userAlphaSpending[1:(kMax - 1)] == 0)) {
message("Changed type of design to ", sQuote(C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY))
typeOfDesign <- C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY
}
}
.assertIsSingleLogical(bindingFutility, "bindingFutility")
.assertIsNumericVector(delayedInformation, "delayedInformation", naAllowed = TRUE)
.assertIsInClosedInterval(delayedInformation, "delayedInformation", lower = 0, upper = NULL, naAllowed = TRUE)
if (designClass == C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL) {
design <- TrialDesignInverseNormal(
kMax = kMax, bindingFutility = bindingFutility,
delayedInformation = delayedInformation
)
} else if (designClass == C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL) {
design <- TrialDesignGroupSequential(
kMax = kMax, bindingFutility = bindingFutility,
delayedInformation = delayedInformation
)
} else {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'designClass' ('", designClass, "') must be '", C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL, "' or ",
"'", C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL, "'"
)
}
.assertIsSingleInteger(sided, "sided", naAllowed = FALSE, validateType = FALSE)
if (!is.integer(sided) && sided %in% c(1, 2)) {
sided <- as.integer(sided)
}
.assertIsSingleCharacter(optimizationCriterion, "optimizationCriterion")
.assertIsSingleCharacter(typeBetaSpending, "typeBetaSpending")
.assertIsSingleLogical(twoSidedPower, "twoSidedPower", naAllowed = TRUE)
.assertIsSingleLogical(betaAdjustment, "betaAdjustment", naAllowed = TRUE)
.assertIsSingleNumber(alpha, "alpha", naAllowed = TRUE)
.assertIsSingleNumber(beta, "beta", naAllowed = TRUE)
.assertIsSingleNumber(deltaWT, "deltaWT", naAllowed = TRUE)
.assertIsSingleNumber(deltaPT1, "deltaPT1", naAllowed = TRUE)
.assertIsSingleNumber(deltaPT0, "deltaPT0", naAllowed = TRUE)
.assertIsSingleNumber(gammaA, "gammaA", naAllowed = TRUE)
.assertIsSingleNumber(gammaB, "gammaB", naAllowed = TRUE)
.assertIsNumericVector(futilityBounds, "futilityBounds", naAllowed = TRUE)
.assertIsNumericVector(informationRates, "informationRates", naAllowed = TRUE)
.assertIsNumericVector(userAlphaSpending, "userAlphaSpending", naAllowed = TRUE)
.assertIsNumericVector(userBetaSpending, "userBetaSpending", naAllowed = TRUE)
design$alpha <- alpha
design$beta <- beta
design$sided <- sided
design$typeOfDesign <- typeOfDesign
design$deltaWT <- deltaWT
design$deltaPT1 <- deltaPT1
design$deltaPT0 <- deltaPT0
design$gammaA <- gammaA
design$gammaB <- gammaB
design$optimizationCriterion <- optimizationCriterion
design$typeBetaSpending <- typeBetaSpending
design$futilityBounds <- futilityBounds
design$informationRates <- informationRates
design$userAlphaSpending <- userAlphaSpending
design$userBetaSpending <- userBetaSpending
design$bindingFutility <- bindingFutility
design$delayedInformation <- delayedInformation
if (!all(is.na(delayedInformation)) && any(delayedInformation > 0)) {
design$.setParameterType("delayedInformation", C_PARAM_USER_DEFINED)
}
if (design$typeOfDesign != C_TYPE_OF_DESIGN_WT_OPTIMUM && optimizationCriterion != C_OPTIMIZATION_CRITERION_DEFAULT) {
warning(
"'optimizationCriterion' (", optimizationCriterion,
") will be ignored because it is only applicable for 'typeOfDesign' = \"", C_TYPE_OF_DESIGN_WT_OPTIMUM, "\"",
call. = FALSE
)
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_HP) {
.assertIsSingleNumber(constantBoundsHP, "constantBoundsHP")
.assertIsInClosedInterval(constantBoundsHP, "constantBoundsHP", lower = 2, upper = NULL)
design$constantBoundsHP <- constantBoundsHP
} else if (constantBoundsHP != C_CONST_BOUND_HP_DEFAULT) {
warning(
"'constantBoundsHP' (", constantBoundsHP,
") will be ignored because it is only applicable for 'typeOfDesign' = \"", C_TYPE_OF_DESIGN_HP, "\"",
call. = FALSE
)
}
if (is.na(twoSidedPower)) {
design$twoSidedPower <- C_TWO_SIDED_POWER_DEFAULT
design$.setParameterType("twoSidedPower", C_PARAM_DEFAULT_VALUE)
} else {
design$twoSidedPower <- twoSidedPower
design$.setParameterType("twoSidedPower", ifelse(
twoSidedPower == C_TWO_SIDED_POWER_DEFAULT, C_PARAM_DEFAULT_VALUE, C_PARAM_USER_DEFINED
))
}
if (design$sided == 2 && grepl("^bs", design$typeBetaSpending)) {
if (is.na(betaAdjustment)) {
design$betaAdjustment <- TRUE
design$.setParameterType("betaAdjustment", C_PARAM_DEFAULT_VALUE)
} else {
design$betaAdjustment <- betaAdjustment
design$.setParameterType("betaAdjustment", ifelse(betaAdjustment, C_PARAM_DEFAULT_VALUE, C_PARAM_USER_DEFINED))
}
} else if (!is.na(betaAdjustment)) {
warning(
"'betaAdjustment' (", betaAdjustment,
") will be ignored because it is only applicable for two-sided beta-spending designs",
call. = FALSE
)
}
design$tolerance <- tolerance
return(design)
}
.getDesignGroupSequentialKMax1 <- function(design) {
design$criticalValues <- .getOneMinusQNorm(design$alpha / design$sided)
design$alphaSpent[1] <- design$alpha
return(invisible(design))
}
#
# Wang and Tsiatis design
#
.getDesignGroupSequentialWangAndTsiatis <- function(design) {
if (design$typeOfDesign == C_TYPE_OF_DESIGN_P) {
design$criticalValues <- getDesignGroupSequentialPocockCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance
)
} else if (design$typeOfDesign == C_TYPE_OF_DESIGN_OF) {
design$criticalValues <- getDesignGroupSequentialOBrienAndFlemingCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance
)
} else {
design$criticalValues <- getDesignGroupSequentialDeltaWTCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance,
design$deltaWT
)
}
.calculateAlphaSpent(design)
return(invisible(design))
}
.getDesignGroupSequentialPampallonaTsiatis <- function(design) {
cppResult <- getDesignGroupSequentialPampallonaTsiatisCpp(
design$tolerance, design$beta, design$alpha, design$kMax,
design$deltaPT0, design$deltaPT1, design$informationRates,
design$sided, design$bindingFutility
)
futilityBounds <- cppResult$futilityBounds
criticalValues <- cppResult$criticalValues
probs <- cppResult$probs
if (design$sided == 1) {
design$betaSpent <- cumsum(probs[1, ])
design$power <- cumsum(probs[3, ] - probs[2, ])
} else {
design$betaSpent <- cumsum(probs[3, ] - probs[2, ])
if (design$twoSidedPower) {
design$power <- cumsum(probs[5, ] - probs[4, ] + probs[1, ])
} else {
design$power <- cumsum(probs[5, ] - probs[4, ])
}
}
design$.setParameterType("betaSpent", C_PARAM_GENERATED)
design$.setParameterType("power", C_PARAM_GENERATED)
design$futilityBounds <- futilityBounds[1:(design$kMax - 1)]
design$criticalValues <- criticalValues
design$.setParameterType("futilityBounds", C_PARAM_GENERATED)
design$.setParameterType("criticalValues", C_PARAM_GENERATED)
.calculateAlphaSpent(design)
design$futilityBounds[design$futilityBounds == 0] <- NA_real_
.assertIsValidBetaSpent(design)
return(invisible(design))
}
.calculateAlphaSpent <- function(design) {
if (design$sided == 2) {
if (design$bindingFutility) {
futilityBounds <- design$futilityBounds
futilityBounds[is.na(futilityBounds)] <- 0
decisionMatrix <- matrix(c(
-design$criticalValues, -futilityBounds, 0,
futilityBounds, 0, design$criticalValues
), nrow = 4, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(-design$criticalValues, design$criticalValues), nrow = 2, byrow = TRUE)
}
} else {
if (design$bindingFutility) {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues
), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
rep(C_FUTILITY_BOUNDS_DEFAULT, design$kMax),
design$criticalValues
), nrow = 2, byrow = TRUE)
}
}
tryCatch(
{
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
if (design$sided == 1) {
design$alphaSpent <- cumsum(probs[3, ] - probs[2, ])
} else if (nrow(decisionMatrix) == 2) {
design$alphaSpent <- cumsum(probs[3, ] - probs[2, ] + probs[1, ])
} else {
design$alphaSpent <- cumsum(probs[5, ] - probs[4, ] + probs[1, ])
}
if (!is.na(design$alphaSpent[design$kMax])) {
design$alphaSpent[design$kMax] <- floor(design$alphaSpent[design$kMax] * 1e8) / 1e8
}
design$.setParameterType("alphaSpent", C_PARAM_GENERATED)
},
error = function(e) {
warning("Failed to calculate 'alphaSpent': ", e, call. = FALSE)
}
)
}
#'
#' Haybittle & Peto design
#'
#' @noRd
#'
.getDesignGroupSequentialHaybittleAndPeto <- function(design) {
scale <- .getOneDimensionalRoot(
function(scale) {
design$criticalValues <- c(rep(design$constantBoundsHP, design$kMax - 1), scale)
if (design$sided == 2) {
decisionMatrix <- matrix(c(-design$criticalValues, design$criticalValues), nrow = 2, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ] + probs[1, ]) - design$alpha)
} else {
if (design$bindingFutility) {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues
), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
rep(C_FUTILITY_BOUNDS_DEFAULT, design$kMax),
design$criticalValues
), nrow = 2, byrow = TRUE)
}
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ]) - design$alpha)
}
},
lower = 0, upper = 8, tolerance = design$tolerance,
callingFunctionInformation = ".getDesignGroupSequentialHaybittleAndPeto"
)
design$criticalValues <- c(rep(design$constantBoundsHP, design$kMax - 1), scale)
.calculateAlphaSpent(design)
return(invisible(design))
}
.getOptimumDesign <- function(deltaWT, design) {
scale <- .getOneDimensionalRoot(
function(scale) {
criticalValues <- scale * design$informationRates^(deltaWT - 0.5)
if (design$sided == 2) {
decisionMatrix <- (matrix(c(-criticalValues, criticalValues), nrow = 2, byrow = TRUE))
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ] + probs[1, ]) - design$alpha)
} else {
if (design$bindingFutility) {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
criticalValues
), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
rep(C_FUTILITY_BOUNDS_DEFAULT, design$kMax),
criticalValues
), nrow = 2, byrow = TRUE)
}
probs <- .getGroupSequentialProbabilities(decisionMatrix, design$informationRates)
return(sum(probs[3, ] - probs[2, ]) - design$alpha)
}
},
lower = 0, upper = 5, tolerance = design$tolerance, callingFunctionInformation = ".getOptimumDesign"
)
design$criticalValues <- scale * design$informationRates^(deltaWT - 0.5)
designCharacteristics <- .getDesignCharacteristics(design = design)
y <- NA_integer_
if (design$optimizationCriterion == C_OPTIMIZATION_CRITERION_ASNH1) {
y <- designCharacteristics$averageSampleNumber1
}
if (design$optimizationCriterion == C_OPTIMIZATION_CRITERION_ASNIFH1) {
y <- designCharacteristics$inflationFactor + designCharacteristics$averageSampleNumber1
}
if (design$optimizationCriterion == C_OPTIMIZATION_CRITERION_ASN_SUM) {
y <- designCharacteristics$averageSampleNumber0 +
designCharacteristics$averageSampleNumber01 + designCharacteristics$averageSampleNumber1
}
return(y)
}
#'
#' Optimum design within Wang and Tsiatis class
#'
#' @noRd
#'
.getDesignGroupSequentialWangAndTsiatisOptimum <- function(design) {
.assertDesignParameterExists(design, "optimizationCriterion", C_OPTIMIZATION_CRITERION_DEFAULT)
.assertIsOptimizationCriterion(design$optimizationCriterion)
optimumDesign <- stats::optimize(.getOptimumDesign,
design = design,
interval = c(0, 1), tol = 0.0001
)
design$deltaWT <- round(optimumDesign$minimum, 3)
design$.setParameterType("deltaWT", C_PARAM_GENERATED)
# Recalculation of design characteristics with rounded design$deltaWT
design$criticalValues <- getDesignGroupSequentialDeltaWTCpp(
design$kMax,
design$alpha,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance,
design$deltaWT
)
designCharacteristics <- .getDesignCharacteristics(design = design)
design$power <- designCharacteristics$power
design$.setParameterType("power", C_PARAM_GENERATED)
.calculateAlphaSpent(design)
return(invisible(design))
}
#'
#' Alpha spending approaches
#'
#' @noRd
#'
.getDesignGroupSequentialAlphaSpending <- function(design, userFunctionCallEnabled) {
design$criticalValues <- getDesignGroupSequentialAlphaSpendingCpp(
design$kMax,
design$alpha,
design$gammaA,
design$typeOfDesign,
design$sided,
design$informationRates,
design$bindingFutility,
design$futilityBounds,
design$tolerance
)
.calculateAlphaSpent(design)
return(.getDesignGroupSequentialBetaSpendingApproaches(design, userFunctionCallEnabled))
}
#'
#' User defined alpha spending approach
#'
#' @noRd
#'
.getDesignGroupSequentialUserDefinedAlphaSpending <- function(design, userFunctionCallEnabled) {
design$criticalValues <- rep(NA_real_, design$kMax)
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY) {
design$userAlphaSpending <- rep(0, design$kMax)
design$userAlphaSpending[design$kMax] <- design$alpha
}
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY && !design$bindingFutility) {
design$criticalValues[1:(design$kMax - 1)] <- C_QNORM_THRESHOLD
design$criticalValues[design$kMax] <- .getOneMinusQNorm(design$alpha / design$sided)
} else {
design$criticalValues <- getDesignGroupSequentialUserDefinedAlphaSpendingCpp(
design$kMax, design$userAlphaSpending, design$sided,
design$informationRates, design$bindingFutility, design$futilityBounds, design$tolerance
)
if (design$typeOfDesign == C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY) {
design$criticalValues[1:(design$kMax - 1)] <- C_QNORM_THRESHOLD
}
}
.calculateAlphaSpent(design)
return(invisible(.getDesignGroupSequentialBetaSpendingApproaches(design, userFunctionCallEnabled)))
}
#'
#' Only for alpha spending approaches
#'
#' @noRd
#'
.getDesignGroupSequentialBetaSpendingApproaches <- function(design, userFunctionCallEnabled) {
# beta spending approaches (additional to alpha spending)!
if (.isBetaSpendingDesignType(design$typeBetaSpending,
userDefinedBetaSpendingIncluded = FALSE, noneIncluded = FALSE
)) {
.getDesignGroupSequentialBetaSpending(design, userFunctionCallEnabled)
}
# User defined beta spending
if (design$typeBetaSpending == C_TYPE_OF_DESIGN_BS_USER) {
.getDesignGroupSequentialUserDefinedBetaSpending(design)
}
return(invisible(design))
}
.fillVec <- function(vec, n) {
return(c(vec, rep(NA_real_, n - length(vec))))
}
#'
#' Beta spending approaches (additional to alpha spending)
#' Find shift with beta spending such that last critical values coincide
#'
#' @noRd
#'
.getDesignGroupSequentialBetaSpending <- function(design, userFunctionCallEnabled) {
cppResult <- getDesignGroupSequentialBetaSpendingCpp(
design$criticalValues,
design$kMax,
design$userAlphaSpending,
design$userBetaSpending,
design$informationRates,
design$bindingFutility,
design$tolerance,
design$typeOfDesign,
design$typeBetaSpending,
design$gammaA,
design$gammaB,
design$alpha,
design$beta,
design$sided,
design$betaAdjustment,
design$twoSidedPower
)
design$futilityBounds <- cppResult$futilityBounds
design$criticalValues <- cppResult$criticalValues
design$betaSpent <- cppResult$betaSpent
design$power <- cppResult$power
if (design$sided == 2) {
.calculateAlphaSpent(design)
}
design$.setParameterType("betaSpent", C_PARAM_GENERATED)
design$.setParameterType("power", C_PARAM_GENERATED)
design$.setParameterType("futilityBounds", C_PARAM_GENERATED)
return(invisible(design))
}
#'
#' User defined beta spending
#'
#' Find shift with beta spending such that last critical values coincide
#'
#' @noRd
#'
.getDesignGroupSequentialUserDefinedBetaSpending <- function(design) {
if (design$typeBetaSpending != C_TYPE_OF_DESIGN_BS_USER) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'typeBetaSpending' ('", design$typeBetaSpending, "') must be '", C_TYPE_OF_DESIGN_BS_USER, "'"
)
}
cppResult <- getDesignGroupSequentialUserDefinedBetaSpendingCpp(
design$criticalValues,
design$kMax,
design$userAlphaSpending,
design$userBetaSpending,
design$sided,
design$informationRates,
design$bindingFutility,
design$tolerance,
design$typeOfDesign,
design$gammaA,
design$alpha,
design$betaAdjustment,
design$twoSidedPower
)
design$futilityBounds <- cppResult$futilityBounds
design$criticalValues <- cppResult$criticalValues
design$betaSpent <- cppResult$betaSpent
design$power <- cppResult$power
if (design$sided == 2) {
.calculateAlphaSpent(design)
}
design$.setParameterType("betaSpent", C_PARAM_GENERATED)
design$.setParameterType("power", C_PARAM_GENERATED)
design$.setParameterType("futilityBounds", C_PARAM_GENERATED)
return(invisible(design))
}
#'
#' Calculate stopping, rejection and futility probabilities for delayed response design
#'
#' @noRd
#'
.calculateDecisionProbabilities <- function(sqrtShift, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues) {
kMax <- length(informationRates)
power <- numeric(kMax)
futilityProbabilities <- numeric(kMax - 1)
stoppingProbabilities <- numeric(kMax - 1)
rejectionProbabilities <- numeric(kMax)
contRegionLower <- c(contRegionLower, decisionCriticalValues[kMax])
for (stage in 1:(kMax)) {
if (!is.na(delayedInformation[stage]) && delayedInformation[stage] > 0) {
# information rate vector in case of recruitment stop at 'stage'
informationRatesUponDelay <- c(
informationRates[1:stage],
informationRates[stage] + delayedInformation[stage]
)
if (stage == 1) {
probs1 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionUpper[stage] - sqrtShift * sqrt(informationRatesUponDelay[1]),
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[2]),
C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
-C_UPPER_BOUNDS_DEFAULT, decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[2]),
contRegionLower[stage] - sqrtShift * sqrt(informationRatesUponDelay[1]), C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
rejectionProbabilities[stage] <- probs1[2, stage + 1] - probs1[1, stage + 1] +
probs2[2, stage + 1] - probs2[1, stage + 1]
power[stage] <- rejectionProbabilities[stage]
futilityProbabilities[stage] <- probs2[2, stage]
stoppingProbabilities[stage] <- probs2[2, stage] + 1 - probs1[1, stage]
} else if (stage < kMax) {
probs1 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift * sqrt(informationRatesUponDelay[1:(stage - 1)]),
contRegionUpper[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage]),
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage + 1]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRatesUponDelay[1:(stage - 1)]),
C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift *
sqrt(informationRatesUponDelay[1:(stage - 1)]), -C_UPPER_BOUNDS_DEFAULT,
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage + 1]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRatesUponDelay[1:(stage - 1)]),
contRegionLower[stage] - sqrtShift * sqrt(informationRatesUponDelay[stage]), C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
rejectionProbabilities[stage] <- probs1[2, stage + 1] - probs1[1, stage + 1] +
probs2[2, stage + 1] - probs2[1, stage + 1]
power[stage] <- sum(rejectionProbabilities[1:stage])
futilityProbabilities[stage] <- probs2[2, stage]
stoppingProbabilities[stage] <- probs2[2, stage] + probs1[2, stage] - probs1[1, stage]
} else {
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
decisionCriticalValues[stage] - sqrtShift * sqrt(informationRates[stage]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRates
)
rejectionProbabilities[stage] <- probs[2, stage] - probs[1, stage]
power[stage] <- sum(rejectionProbabilities[1:stage])
}
} else {
if (stage == 1) {
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[stage] - sqrtShift * sqrt(informationRates[stage]),
contRegionUpper[stage] - sqrtShift * sqrt(informationRates[stage])
),
nrow = 2, byrow = TRUE
), informationRates[1]
)
} else {
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
contRegionLower[stage] - sqrtShift * sqrt(informationRates[stage]),
contRegionUpper[1:(stage - 1)] - sqrtShift * sqrt(informationRates[1:(stage - 1)]),
contRegionUpper[stage] - sqrtShift * sqrt(informationRates[stage])
),
nrow = 2, byrow = TRUE
), informationRates[1:stage]
)
}
rejectionProbabilities[stage] <- probs[3, stage] - probs[2, stage]
if (stage < kMax) {
futilityProbabilities[stage] <- probs[1, stage]
stoppingProbabilities[stage] <- futilityProbabilities[stage] + rejectionProbabilities[stage]
}
power[stage] <- sum(rejectionProbabilities[1:stage])
}
}
return(list(
probs = probs,
power = power,
rejectionProbabilities = rejectionProbabilities,
futilityProbabilities = futilityProbabilities,
stoppingProbabilities = stoppingProbabilities
))
}
#'
#' @title
#' Get Design Inverse Normal
#'
#' @description
#' Provides adjusted boundaries and defines a group sequential design for its use in
#' the inverse normal combination test.
#'
#' @inheritParams getDesignGroupSequential
#'
#' @template details_group_sequential_design
#'
#' @template return_object_trial_design
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_design_inverse_normal
#'
#' @export
#'
getDesignInverseNormal <- function(...,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = 1L, # C_SIDED_DEFAULT
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = c("OF", "P", "WT", "PT", "HP", "WToptimum", "asP", "asOF", "asKD", "asHSD", "asUser", "noEarlyEfficacy"), # C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = c("ASNH1", "ASNIFH1", "ASNsum"), # C_OPTIMIZATION_CRITERION_DEFAULT
gammaA = NA_real_,
typeBetaSpending = c("none", "bsP", "bsOF", "bsKD", "bsHSD", "bsUser"), # C_TYPE_OF_DESIGN_BS_NONE
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = NA,
betaAdjustment = NA,
constantBoundsHP = 3, # C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
tolerance = 1e-08 # C_DESIGN_TOLERANCE_DEFAULT
) {
.warnInCaseOfUnknownArguments(functionName = "getDesignInverseNormal", ...)
return(.getDesignGroupSequential(
designClass = C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
betaAdjustment = betaAdjustment,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
tolerance = tolerance,
userFunctionCallEnabled = TRUE
))
}
.getDesignInverseNormal <- function(...,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = NA,
betaAdjustment = NA,
constantBoundsHP = C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
tolerance = C_DESIGN_TOLERANCE_DEFAULT) {
.warnInCaseOfUnknownArguments(functionName = "getDesignInverseNormal", ...)
return(.getDesignGroupSequential(
designClass = C_CLASS_NAME_TRIAL_DESIGN_INVERSE_NORMAL,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
betaAdjustment = betaAdjustment,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
tolerance = tolerance,
userFunctionCallEnabled = FALSE
))
}
.getDesignGroupSequentialDefaultValues <- function() {
return(list(
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
twoSidedPower = C_TWO_SIDED_POWER_DEFAULT,
tolerance = C_DESIGN_TOLERANCE_DEFAULT
))
}
.getDesignInverseNormalDefaultValues <- function() {
return(.getDesignGroupSequentialDefaultValues())
}
#
# Param: userFunctionCallEnabled if \code{TRUE}, additional parameter validation methods will be called.
#
.getDesignGroupSequential <- function(...,
designClass = C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = C_SIDED_DEFAULT,
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = C_OPTIMIZATION_CRITERION_DEFAULT,
gammaA = NA_real_,
typeBetaSpending = C_TYPE_OF_DESIGN_BS_NONE,
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = C_BINDING_FUTILITY_DEFAULT,
constantBoundsHP = C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
betaAdjustment = NA,
delayedInformation = NA_real_,
tolerance = C_DESIGN_TOLERANCE_DEFAULT,
userFunctionCallEnabled = FALSE) {
typeOfDesign <- .matchArgument(typeOfDesign, C_DEFAULT_TYPE_OF_DESIGN)
optimizationCriterion <- .matchArgument(optimizationCriterion, C_OPTIMIZATION_CRITERION_DEFAULT)
typeBetaSpending <- .matchArgument(typeBetaSpending, C_TYPE_OF_DESIGN_BS_NONE)
if (.isDefinedArgument(kMax, argumentExistsValidationEnabled = userFunctionCallEnabled)) {
.assertIsValidKMax(kMax, showWarnings = TRUE)
if (!is.integer(kMax)) {
kMax <- as.integer(kMax)
}
}
if (is.na(bindingFutility)) {
bindingFutility <- C_BINDING_FUTILITY_DEFAULT
} else if (userFunctionCallEnabled && typeOfDesign != C_TYPE_OF_DESIGN_PT &&
!(typeBetaSpending == "bsP" || typeBetaSpending == "bsOF" || typeBetaSpending == "bsKD" ||
typeBetaSpending == "bsHSD" || typeBetaSpending == "bsUser") &&
((!is.na(kMax) && kMax == 1) || any(is.na(futilityBounds)) ||
(!any(is.na(futilityBounds)) && all(futilityBounds == C_FUTILITY_BOUNDS_DEFAULT)))) {
warning("'bindingFutility' (", bindingFutility, ") will be ignored", call. = FALSE)
}
design <- .createDesign(
designClass = designClass,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
betaAdjustment = betaAdjustment,
delayedInformation = delayedInformation,
tolerance = tolerance
)
if (userFunctionCallEnabled) {
.validateBaseParameters(design, twoSidedWarningForDefaultValues = FALSE)
.validateTypeOfDesign(design)
.assertIsValidTolerance(tolerance)
.assertDesignParameterExists(design, "alpha", C_ALPHA_DEFAULT)
.assertDesignParameterExists(design, "beta", C_BETA_DEFAULT)
.assertDesignParameterExists(design, "sided", C_SIDED_DEFAULT)
.assertDesignParameterExists(design, "typeOfDesign", C_DEFAULT_TYPE_OF_DESIGN)
.assertDesignParameterExists(design, "bindingFutility", C_BINDING_FUTILITY_DEFAULT)
.assertDesignParameterExists(design, "tolerance", C_DESIGN_TOLERANCE_DEFAULT)
if (typeOfDesign != C_TYPE_OF_DESIGN_PT) {
if (!is.na(deltaPT1)) {
warning("'deltaPT1' (", deltaPT1, ") will be ignored", call. = FALSE)
}
if (!is.na(deltaPT0)) {
warning("'deltaPT0' (", deltaPT0, ") will be ignored", call. = FALSE)
}
}
if (typeOfDesign != C_TYPE_OF_DESIGN_WT && !is.na(deltaWT)) {
warning("'deltaWT' (", deltaWT, ") will be ignored", call. = FALSE)
}
if (typeOfDesign != C_TYPE_OF_DESIGN_AS_KD &&
typeOfDesign != C_TYPE_OF_DESIGN_AS_HSD && !is.na(gammaA)) {
warning("'gammaA' (", gammaA, ") will be ignored", call. = FALSE)
}
if (typeBetaSpending != C_TYPE_OF_DESIGN_BS_KD &&
typeBetaSpending != C_TYPE_OF_DESIGN_BS_HSD && !is.na(gammaB)) {
warning("'gammaB' (", gammaB, ") will be ignored", call. = FALSE)
}
if (typeBetaSpending != C_TYPE_OF_DESIGN_BS_USER && !all(is.na(userBetaSpending))) {
warning("'userBetaSpending' (", .arrayToString(userBetaSpending), ") will be ignored", call. = FALSE)
}
if (!(typeOfDesign %in% c(C_TYPE_OF_DESIGN_AS_USER)) &&
!all(is.na(userAlphaSpending))) {
warning("'userAlphaSpending' (", .arrayToString(userAlphaSpending), ") will be ignored", call. = FALSE)
}
}
if (design$sided == 2 && design$bindingFutility && design$typeOfDesign != C_TYPE_OF_DESIGN_PT &&
!.isBetaSpendingDesignType(design$typeBetaSpending)) {
warning("'bindingFutility' will be ignored because the test is defined as two-sided", call. = FALSE)
design$bindingFutility <- FALSE
}
if (design$sided == 1 && design$twoSidedPower) {
warning("'twoSidedPower' will be ignored because the test is defined as one-sided", call. = FALSE)
design$twoSidedPower <- FALSE
}
if (userFunctionCallEnabled) {
.validateAlphaAndBeta(design)
}
design$alphaSpent <- rep(NA_real_, design$kMax)
design$betaSpent <- rep(NA_real_, design$kMax)
design$power <- rep(NA_real_, design$kMax)
if (userFunctionCallEnabled) {
design$.setParameterType("betaSpent", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("power", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("alphaSpent", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("criticalValues", C_PARAM_GENERATED)
}
if (design$kMax == 1) {
.getDesignGroupSequentialKMax1(design)
} else {
# Wang and Tsiatis design
if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT ||
design$typeOfDesign == C_TYPE_OF_DESIGN_P ||
design$typeOfDesign == C_TYPE_OF_DESIGN_OF) {
.getDesignGroupSequentialWangAndTsiatis(design)
}
# Pampallona & Tsiatis design
else if (design$typeOfDesign == C_TYPE_OF_DESIGN_PT) {
.getDesignGroupSequentialPampallonaTsiatis(design)
}
# Haybittle & Peto design
else if (design$typeOfDesign == C_TYPE_OF_DESIGN_HP) {
.getDesignGroupSequentialHaybittleAndPeto(design)
}
# Optimum design within Wang and Tsiatis class
else if (design$typeOfDesign == C_TYPE_OF_DESIGN_WT_OPTIMUM) {
.getDesignGroupSequentialWangAndTsiatisOptimum(design)
}
# alpha spending approaches
else if (.isAlphaSpendingDesignType(design$typeOfDesign, userDefinedAlphaSpendingIncluded = FALSE)) {
.getDesignGroupSequentialAlphaSpending(design, userFunctionCallEnabled)
}
# user defined alpha spending approach
else if (design$typeOfDesign %in% c(C_TYPE_OF_DESIGN_AS_USER, C_TYPE_OF_DESIGN_NO_EARLY_EFFICACY)) {
.getDesignGroupSequentialUserDefinedAlphaSpending(design, userFunctionCallEnabled)
} else {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "no calculation routine defined for ", design$typeOfDesign)
}
}
design$stageLevels <- 1 - stats::pnorm(design$criticalValues)
design$.setParameterType("stageLevels", C_PARAM_GENERATED)
if (design$kMax == 1) {
design$.setParameterType("futilityBounds", C_PARAM_NOT_APPLICABLE)
}
if (!all(is.na(design$futilityBounds))) {
if (length(design$futilityBounds) == 0 || all(design$futilityBounds == C_FUTILITY_BOUNDS_DEFAULT)) {
design$.setParameterType("bindingFutility", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("futilityBounds", C_PARAM_NOT_APPLICABLE)
} else if (userFunctionCallEnabled &&
any(design$futilityBounds > design$criticalValues[1:(design$kMax - 1)] - 0.01, na.rm = TRUE)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'futilityBounds' (", .arrayToString(design$futilityBounds), ") too extreme for this situation"
)
}
}
.assertIsValidAlphaSpent(design, userFunctionCallEnabled)
design$.initStages()
# we use 7.5 instead of C_QNORM_THRESHOLD as threshold
design$criticalValues[!is.na(design$criticalValues) & design$criticalValues <= -7.5] <- -Inf
design$criticalValues[!is.na(design$criticalValues) & design$criticalValues >= 7.5] <- Inf
design$futilityBounds[!is.na(design$futilityBounds) & design$futilityBounds <=
C_FUTILITY_BOUNDS_DEFAULT] <- C_FUTILITY_BOUNDS_DEFAULT
if (design$kMax == 1) {
if (!identical(design$informationRates, 1)) {
if (!is.na(design$informationRates)) {
warning("Information rate", ifelse(length(design$informationRates) != 1, "s", ""), " ",
.arrayToString(design$informationRates, vectorLookAndFeelEnabled = TRUE),
" will be ignored",
call. = FALSE
)
}
design$informationRates <- 1
}
design$.setParameterType("informationRates", C_PARAM_NOT_APPLICABLE)
design$.setParameterType("stages", C_PARAM_NOT_APPLICABLE)
}
.assertIsNumericVector(delayedInformation, "delayedInformation", naAllowed = TRUE)
if (all(is.na(delayedInformation))) {
# delayed response design is disabled
return(design)
}
if (all(!is.na(delayedInformation)) && all(delayedInformation < 1e-03)) {
warning("At least one delayed information value must be >= 1e-03 to enable delayed response.",
" 'delayedInformation' (", .arrayToString(delayedInformation), ") will be ignored",
call. = FALSE
)
return(design)
}
# proceed with delayed response design
.assertIsInClosedInterval(delayedInformation, "delayedInformation", lower = 0, upper = NULL)
kMax <- design$kMax
contRegionUpper <- design$criticalValues
contRegionLower <- design$futilityBounds
informationRates <- design$informationRates
decisionCriticalValues <- numeric(kMax)
reversalProbabilities <- numeric(kMax - 1)
if (!all(is.na(delayedInformation)) &&
(design$sided != 1 || all(design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT + 1e-06))) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"decision critical values for delayed response design are only available for ",
"one-sided designs with valid futility bounds"
)
}
if (length(delayedInformation) == 1) {
delayedInformation <- rep(delayedInformation, kMax - 1)
}
if (length(delayedInformation) != kMax - 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'delayedInformation' (", .arrayToString(delayedInformation), ") must have length ",
(kMax - 1), " (kMax - 1)"
)
}
indices <- which(delayedInformation > 0 & delayedInformation < 1e-03)
n <- length(indices)
if (n > 0) {
warning("The", ifelse(n == 1, "", paste0(" ", n)), " delayed information value", ifelse(n == 1, "", "s"), " ",
.arrayToString(delayedInformation[indices], mode = "and"),
" will be replaced by 1e-03 to achieve reasonable results",
call. = FALSE
)
delayedInformation[indices] <- 1e-03
}
# sensible interim choices are restricted by amount of delayed information
eps <- design$informationRates[1:(design$kMax - 1)] + delayedInformation
if (!any(is.na(eps)) && any(eps >= 1)) {
stages <- which(eps >= 1)
stagesInfo <- ifelse(length(stages) == 1, paste0(" ", stages), paste0("s ", .arrayToString(stages, mode = "and")))
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'delayedInformation[stage] + informationRates[stage]' for ",
"stage", stagesInfo, " too large (>= 1). Recruitment stop analysis information + pipeline data ",
"information cannot exceed overall trial information. Instead, the ",
"recruitment stop analysis would be skipped, directly proceeding to the ",
"final analysis"
)
}
# loop iterating through the stages calculation the decision critical values
for (stage in 1:(kMax - 1)) {
if (!is.na(delayedInformation[stage]) && delayedInformation[stage] >= 1e-03 - 1e-06) {
# information rate vector in case of recruitment stop at 'stage'
informationRatesUponDelay <- c(
informationRates[1:stage],
informationRates[stage] + delayedInformation[stage]
)
if (stage == 1) {
decisionCriticalValues[stage] <- .getOneDimensionalRoot(function(secondCriticalValue) {
probs1 <- .getGroupSequentialProbabilities(
matrix(c(contRegionUpper[stage], secondCriticalValue, C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
-C_UPPER_BOUNDS_DEFAULT, secondCriticalValue,
contRegionLower[stage], C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
return(probs1[1, stage + 1] - probs2[2, stage + 1] + probs2[1, stage + 1])
}, lower = -C_UPPER_BOUNDS_DEFAULT, upper = C_UPPER_BOUNDS_DEFAULT, tolerance = design$tolerance)
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionUpper[stage], decisionCriticalValues[stage],
C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
} else {
decisionCriticalValues[stage] <- .getOneDimensionalRoot(function(secondCriticalValue) {
probs1 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)], contRegionUpper[stage], secondCriticalValue,
contRegionUpper[1:(stage - 1)], C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
probs2 <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)], -C_UPPER_BOUNDS_DEFAULT, secondCriticalValue,
contRegionUpper[1:(stage - 1)], contRegionLower[stage], C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
return(probs1[1, stage + 1] - probs2[2, stage + 1] + probs2[1, stage + 1])
}, lower = -C_UPPER_BOUNDS_DEFAULT, upper = C_UPPER_BOUNDS_DEFAULT, tolerance = design$tolerance)
probs <- .getGroupSequentialProbabilities(
matrix(
c(
contRegionLower[1:(stage - 1)], contRegionUpper[stage], decisionCriticalValues[stage],
contRegionUpper[1:(stage - 1)], C_UPPER_BOUNDS_DEFAULT, C_UPPER_BOUNDS_DEFAULT
),
nrow = 2, byrow = TRUE
), informationRatesUponDelay
)
}
if (stage < kMax) {
reversalProbabilities[stage] <- probs[1, stage + 1]
}
} else {
decisionCriticalValues[stage] <- NA_real_
reversalProbabilities[stage] <- NA_real_
}
decisionCriticalValues[kMax] <- contRegionUpper[kMax]
alphaSpent <- .calculateDecisionProbabilities(
sqrtShift = 0, informationRates,
delayedInformation, contRegionUpper, contRegionLower, decisionCriticalValues
)$power
}
decisionCriticalValues[decisionCriticalValues <= -C_UPPER_BOUNDS_DEFAULT + 1e-06] <- NA_real_
decisionCriticalValues[decisionCriticalValues >= C_UPPER_BOUNDS_DEFAULT - 1e-06] <- NA_real_
design$decisionCriticalValues <- decisionCriticalValues
design$reversalProbabilities <- reversalProbabilities
design$delayedInformation <- delayedInformation
design$delayedInformation[design$delayedInformation < 1e-03] <- 0
design$.setParameterType("decisionCriticalValues", C_PARAM_GENERATED)
design$.setParameterType("reversalProbabilities", C_PARAM_GENERATED)
warning("The delayed information design feature is experimental and ",
"hence not fully validated (see www.rpact.com/experimental)",
call. = FALSE
)
return(design)
}
# to avoid error messages in case of repeated p-values computation
.assertIsValidAlphaSpent <- function(design, userFunctionCallEnabled = TRUE) {
if (!userFunctionCallEnabled) {
return(invisible())
}
if (design$informationRates[design$kMax] != 1) {
return(invisible())
}
if (is.na(design$alphaSpent[design$kMax]) || abs(design$alphaSpent[design$kMax] - design$alpha) > 1e-05) {
stop(
C_EXCEPTION_TYPE_RUNTIME_ISSUE, "critical values cannot be calculated ",
"(alpha: ", design$alpha, "; alpha spent at maximum stage: ", design$alphaSpent[design$kMax], ")"
)
}
}
.assertIsValidBetaSpent <- function(design, ..., userFunctionCallEnabled = TRUE, iteration = 1) {
if (!userFunctionCallEnabled) {
return(invisible())
}
if (design$informationRates[design$kMax] != 1) {
return(invisible())
}
if (iteration < 0) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "critical values cannot be calculated")
}
if (is.na(design$betaSpent[design$kMax]) || abs(design$betaSpent[design$kMax] - design$beta) > 1e-05) {
stop(
C_EXCEPTION_TYPE_RUNTIME_ISSUE, "critical values cannot be calculated ",
"(beta spent at maximum stage: ", design$betaSpent[design$kMax], ")"
)
}
}
#'
#' @title
#' Get Design Group Sequential
#'
#' @description
#' Provides adjusted boundaries and defines a group sequential design.
#'
#' @inheritParams param_kMax
#' @inheritParams param_alpha
#' @inheritParams param_beta
#' @inheritParams param_sided
#' @inheritParams param_typeOfDesign
#' @inheritParams param_informationRates
#' @param futilityBounds The futility bounds, defined on the test statistic z scale
#' (numeric vector of length \code{kMax - 1}).
#' @inheritParams param_bindingFutility
#' @param deltaWT Delta for Wang & Tsiatis Delta class.
#' @param deltaPT1 Delta1 for Pampallona & Tsiatis class rejecting H0 boundaries.
#' @param deltaPT0 Delta0 for Pampallona & Tsiatis class rejecting H1 boundaries.
#' @param constantBoundsHP The constant bounds up to stage \code{kMax - 1} for the
#' Haybittle & Peto design (default is \code{3}).
#' @param optimizationCriterion Optimization criterion for optimum design within
#' Wang & Tsiatis class (\code{"ASNH1"}, \code{"ASNIFH1"},
#' \code{"ASNsum"}), default is \code{"ASNH1"}, see details.
#' @param typeBetaSpending Type of beta spending. Type of of beta spending is one of the following:
#' O'Brien & Fleming type beta spending, Pocock type beta spending,
#' Kim & DeMets beta spending, Hwang, Shi & DeCani beta spending, user defined
#' beta spending (\code{"bsOF"}, \code{"bsP"}, \code{"bsKD"},
#' \code{"bsHSD"}, \code{"bsUser"}, default is \code{"none"}).
#' @param gammaA Parameter for alpha spending function.
#' @param gammaB Parameter for beta spending function.
#' @inheritParams param_userAlphaSpending
#' @param delayedInformation Delay of information for delayed response designs. Can be a numeric value or a
#' numeric vector of length \code{kMax - 1}
#' @param userBetaSpending The user defined beta spending. Vector of length \code{kMax} containing the cumulative
#' beta-spending up to each interim stage.
#' @param twoSidedPower For two-sided testing, if \code{twoSidedPower = TRUE} is specified
#' the sample size calculation is performed by considering both tails of the distribution.
#' Default is \code{FALSE}, i.e., it is assumed that one tail probability is equal to 0 or the power
#' should be directed to one part.
#' @param betaAdjustment For two-sided beta spending designs, if \code{betaAdjustement = TRUE} a linear
#' adjustment of the beta spending values is performed if an overlapping of decision regions for futility
#' stopping at earlier stages occurs, otherwise no adjustment is performed (default is \code{TRUE}).
#' @param tolerance The numerical tolerance, default is \code{1e-08}.
#' @inheritParams param_three_dots
#'
#' @template details_group_sequential_design
#'
#' @template return_object_trial_design
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_design_group_sequential
#'
#' @export
#'
getDesignGroupSequential <- function(...,
kMax = NA_integer_,
alpha = NA_real_,
beta = NA_real_,
sided = 1L, # C_SIDED_DEFAULT
informationRates = NA_real_,
futilityBounds = NA_real_,
typeOfDesign = c("OF", "P", "WT", "PT", "HP", "WToptimum", "asP", "asOF", "asKD", "asHSD", "asUser", "noEarlyEfficacy"), # C_DEFAULT_TYPE_OF_DESIGN,
deltaWT = NA_real_,
deltaPT1 = NA_real_,
deltaPT0 = NA_real_,
optimizationCriterion = c("ASNH1", "ASNIFH1", "ASNsum"), # C_OPTIMIZATION_CRITERION_DEFAULT
gammaA = NA_real_,
typeBetaSpending = c("none", "bsP", "bsOF", "bsKD", "bsHSD", "bsUser"), # C_TYPE_OF_DESIGN_BS_NONE
userAlphaSpending = NA_real_,
userBetaSpending = NA_real_,
gammaB = NA_real_,
bindingFutility = NA,
betaAdjustment = NA,
constantBoundsHP = 3, # C_CONST_BOUND_HP_DEFAULT,
twoSidedPower = NA,
delayedInformation = NA_real_,
tolerance = 1e-08 # C_DESIGN_TOLERANCE_DEFAULT
) {
.warnInCaseOfUnknownArguments(functionName = "getDesignGroupSequential", ...)
return(.getDesignGroupSequential(
designClass = C_CLASS_NAME_TRIAL_DESIGN_GROUP_SEQUENTIAL,
kMax = kMax,
alpha = alpha,
beta = beta,
sided = sided,
informationRates = informationRates,
futilityBounds = futilityBounds,
typeOfDesign = typeOfDesign,
deltaWT = deltaWT,
deltaPT1 = deltaPT1,
deltaPT0 = deltaPT0,
optimizationCriterion = optimizationCriterion,
gammaA = gammaA,
typeBetaSpending = typeBetaSpending,
userAlphaSpending = userAlphaSpending,
userBetaSpending = userBetaSpending,
gammaB = gammaB,
bindingFutility = bindingFutility,
constantBoundsHP = constantBoundsHP,
twoSidedPower = twoSidedPower,
betaAdjustment = betaAdjustment,
delayedInformation = delayedInformation,
tolerance = tolerance,
userFunctionCallEnabled = TRUE
))
}
.getFixedSampleSize <- function(alpha, beta, sided, twoSidedPower = C_TWO_SIDED_POWER_DEFAULT) {
.assertIsValidAlphaAndBeta(alpha = alpha, beta = beta)
.assertIsValidSidedParameter(sided)
if (sided == 1) {
return((.getOneMinusQNorm(alpha) + .getOneMinusQNorm(beta))^2)
}
if (twoSidedPower) {
n <- .getOneDimensionalRoot(
function(n) {
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n)) -
stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n)) + beta
},
lower = 0,
upper = 2 * (.getOneMinusQNorm(alpha / 2) + .getOneMinusQNorm(beta))^2,
tolerance = 1e-08, callingFunctionInformation = ".getFixedSampleSize"
)
} else {
n <- (.getOneMinusQNorm(alpha / 2) + .getOneMinusQNorm(beta))^2
}
return(n)
}
#' @title
#' Get Design Characteristics
#'
#' @description
#' Calculates the characteristics of a design and returns it.
#'
#' @inheritParams param_design
#' @inheritParams param_three_dots
#'
#' @details
#' Calculates the inflation factor (IF),
#' the expected reduction in sample size under H1, under H0, and under a value in between H0 and H1.
#' Furthermore, absolute information values are calculated
#' under the prototype case testing H0: mu = 0 against H1: mu = 1.
#'
#' @return Returns a \code{\link{TrialDesignCharacteristics}} object.
#' The following generics (R generic functions) are available for this result object:
#' \itemize{
#' \item \code{\link[=names.FieldSet]{names()}} to obtain the field names,
#' \item \code{\link[=print.FieldSet]{print()}} to print the object,
#' \item \code{\link[=summary.ParameterSet]{summary()}} to display a summary of the object,
#' \item \code{\link[=plot.ParameterSet]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.TrialDesignCharacteristics]{as.data.frame()}} to coerce the object to a \code{\link[base]{data.frame}},
#' \item \code{\link[=as.matrix.FieldSet]{as.matrix()}} to coerce the object to a \code{\link[base]{matrix}}.
#' }
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_design_characteristics
#'
#' @export
#'
getDesignCharacteristics <- function(design = NULL, ...) {
if (is.null(design)) {
design <- .getDefaultDesign(..., type = "characteristics")
.warnInCaseOfUnknownArguments(
functionName = "getDesignCharacteristics",
ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = FALSE), ...
)
} else {
.assertIsTrialDesign(design)
.warnInCaseOfUnknownArguments(functionName = "getDesignCharacteristics", ...)
.warnInCaseOfTwoSidedPowerArgument(...)
}
return(.getDesignCharacteristics(design = design, userFunctionCallEnabled = TRUE))
}
.getDesignCharacteristics <- function(..., design, userFunctionCallEnabled = FALSE) {
.assertIsTrialDesignInverseNormalOrGroupSequential(design)
.assertDesignParameterExists(design, "sided", C_SIDED_DEFAULT)
.assertIsValidSidedParameter(design$sided)
if (userFunctionCallEnabled) {
.validateAlphaAndBeta(design = design)
}
design$informationRates <- .getValidatedInformationRates(design, writeToDesign = FALSE)
if ((design$typeOfDesign == C_TYPE_OF_DESIGN_PT ||
.isBetaSpendingDesignType(design$typeBetaSpending)) && design$sided == 2 && design$kMax == 2) {
design$futilityBounds[is.na(design$futilityBounds)] <- 0
}
design$futilityBounds <- .getValidatedFutilityBounds(design,
writeToDesign = FALSE, twoSidedWarningForDefaultValues = FALSE
)
designCharacteristics <- TrialDesignCharacteristics(design = design)
designCharacteristics$rejectionProbabilities <- rep(NA_real_, design$kMax)
designCharacteristics$.setParameterType("rejectionProbabilities", C_PARAM_NOT_APPLICABLE)
designCharacteristics$futilityProbabilities <- rep(NA_real_, design$kMax - 1)
designCharacteristics$.setParameterType("futilityProbabilities", C_PARAM_NOT_APPLICABLE)
nFixed <- .getFixedSampleSize(
alpha = design$alpha, beta = design$beta,
sided = design$sided, twoSidedPower = design$twoSidedPower
)
designCharacteristics$nFixed <- nFixed
designCharacteristics$.setParameterType("nFixed", C_PARAM_GENERATED)
design$criticalValues[design$criticalValues > 7.5] <- 7.5
if (length(design$decisionCriticalValues) > 0) {
design$decisionCriticalValues[!is.na(design$decisionCriticalValues) & design$decisionCriticalValues > 7.5] <- 7.5
}
informationRates <- design$informationRates
if (design$kMax == 1) {
designCharacteristics$shift <- nFixed
designCharacteristics$.setParameterType("shift", C_PARAM_GENERATED)
designCharacteristics$inflationFactor <- designCharacteristics$shift / nFixed
designCharacteristics$.setParameterType("inflationFactor", C_PARAM_GENERATED)
designCharacteristics$power <- 1 - design$beta
designCharacteristics$.setParameterType("power", design$.getParameterType("power"))
designCharacteristics$.setParameterType("information", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType("averageSampleNumber1", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType("averageSampleNumber01", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType("averageSampleNumber0", C_PARAM_NOT_APPLICABLE)
designCharacteristics$.setParameterType(".probs", C_PARAM_NOT_APPLICABLE)
return(designCharacteristics)
}
if (!any(is.na(design$delayedInformation)) && length(design$decisionCriticalValues) > 0) {
kMax <- design$kMax
contRegionUpper <- design$criticalValues
contRegionLower <- design$futilityBounds
informationRates <- design$informationRates
decisionCriticalValues <- design$decisionCriticalValues
informationRates <- design$informationRates
delayedInformation <- design$delayedInformation
kMax <- length(informationRates)
shift <- .getOneDimensionalRoot(
function(shift) {
resultsH1 <- .calculateDecisionProbabilities(sqrt(shift), informationRates, delayedInformation, contRegionUpper, contRegionLower, decisionCriticalValues)
return(resultsH1$power[kMax] - 1 + design$beta)
},
lower = 0, upper = 4 * nFixed,
tolerance = design$tolerance, callingFunctionInformation = ".getDesignCharacteristics"
)
stopping <- numeric(kMax)
futility <- numeric(kMax)
rejectionProbabilities <- numeric(kMax)
resultsH1 <- .calculateDecisionProbabilities(
sqrtShift = sqrt(shift), informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues
)
resultsH01 <- .calculateDecisionProbabilities(
sqrtShift = sqrt(shift) / 2, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues
)
resultsH0 <- .calculateDecisionProbabilities(
sqrtShift = 0, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues
)
designCharacteristics$shift <- shift
designCharacteristics$.probs <- resultsH1$probs
designCharacteristics$power <- resultsH1$power
designCharacteristics$information <- informationRates * shift
designCharacteristics$averageSampleNumber1 <-
(shift - sum(resultsH1$stoppingProbabilities * (informationRates[kMax] -
delayedInformation - informationRates[1:(kMax - 1)]) * shift)) / nFixed
designCharacteristics$averageSampleNumber01 <-
(shift - sum(resultsH01$stoppingProbabilities * (informationRates[kMax] -
delayedInformation - informationRates[1:(kMax - 1)]) * shift)) / nFixed
designCharacteristics$averageSampleNumber0 <-
(shift - sum(resultsH0$stoppingProbabilities * (informationRates[kMax] -
delayedInformation - informationRates[1:(kMax - 1)]) * shift)) / nFixed
futilityProbabilities <- resultsH1$futilityProbabilities
rejectionProbabilities <- resultsH1$power
stoppingProbabilities <- resultsH1$stoppingProbabilities
rejectionProbabilities[2:kMax] <- resultsH1$power[2:kMax] - rejectionProbabilities[1:(kMax - 1)]
designCharacteristics$rejectionProbabilities <- rejectionProbabilities
designCharacteristics$futilityProbabilities <- futilityProbabilities
} else if ((design$typeOfDesign == C_TYPE_OF_DESIGN_PT ||
.isBetaSpendingDesignType(design$typeBetaSpending)) && design$sided == 2) {
design$futilityBounds[is.na(design$futilityBounds)] <- 0
shift <- .getOneDimensionalRoot(
function(shift) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
c(-design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
c(design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 4, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
if (design$twoSidedPower) {
return(sum(probs[5, ] - probs[4, ] + probs[1, ]) - 1 + design$beta)
} else {
return(sum(probs[5, ] - probs[4, ]) - 1 + design$beta)
}
},
lower = 0, upper = 4 * (.getOneMinusQNorm(design$alpha / design$sided) + .getOneMinusQNorm(design$beta))^2,
tolerance = design$tolerance, callingFunctionInformation = ".getDesignCharacteristics"
)
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
c(-design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
c(design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]), 0),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 4, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$shift <- shift
designCharacteristics$.probs <- probs
if (design$twoSidedPower) {
designCharacteristics$power <- cumsum(probs[5, ] - probs[4, ] + probs[1, ])
designCharacteristics$rejectionProbabilities <- probs[5, ] - probs[4, ] + probs[1, ]
} else {
designCharacteristics$power <- cumsum(probs[5, ] - probs[4, ])
designCharacteristics$rejectionProbabilities <- probs[5, ] - probs[4, ]
}
if (design$kMax > 1) {
designCharacteristics$futilityProbabilities <- probs[3, 1:(design$kMax - 1)] - probs[2, 1:(design$kMax - 1)]
}
designCharacteristics$information <- informationRates * shift
designCharacteristics$averageSampleNumber1 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs, shift, nFixed
)
decisionMatrix <- matrix(c(
-design$criticalValues,
c(-design$futilityBounds, 0),
c(design$futilityBounds, 0),
design$criticalValues
), nrow = 4, byrow = TRUE)
probs0 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber0 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs0, shift, nFixed
)
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates) / 2,
c(-design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]) / 2, 0),
c(design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]) / 2, 0),
design$criticalValues - sqrt(shift * informationRates) / 2
), nrow = 4, byrow = TRUE)
probs01 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber01 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs01, shift, nFixed
)
design$futilityBounds[design$futilityBounds == 0] <- NA_real_
} else {
shift <- .getOneDimensionalRoot(
function(shift) {
if (design$sided == 2) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
if (design$twoSidedPower) {
return(sum(probs[3, ] - probs[2, ] + probs[1, ]) - 1 + design$beta)
} else {
return(sum(probs[3, ] - probs[2, ]) - 1 + design$beta)
}
} else {
shiftedFutilityBounds <- design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)])
shiftedFutilityBounds[design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <-
C_FUTILITY_BOUNDS_DEFAULT
decisionMatrix <- matrix(c(
shiftedFutilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
return(sum(probs[3, ] - probs[2, ]) - 1 + design$beta)
}
},
lower = 0, upper = 4 * (.getOneMinusQNorm(design$alpha / design$sided) + .getOneMinusQNorm(design$beta))^2,
tolerance = design$tolerance, callingFunctionInformation = ".getDesignCharacteristics"
)
if (design$sided == 2) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates),
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
} else {
shiftedFutilityBounds <- design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)])
shiftedFutilityBounds[design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <-
C_FUTILITY_BOUNDS_DEFAULT
decisionMatrix <- matrix(c(
shiftedFutilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues - sqrt(shift * informationRates)
), nrow = 2, byrow = TRUE)
}
probs <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$shift <- shift
designCharacteristics$.probs <- probs
if (design$twoSidedPower) {
designCharacteristics$power <- cumsum(probs[3, ] - probs[2, ] + probs[1, ])
designCharacteristics$rejectionProbabilities <- probs[3, ] - probs[2, ] + probs[1, ]
} else {
designCharacteristics$power <- cumsum(probs[3, ] - probs[2, ])
designCharacteristics$rejectionProbabilities <- probs[3, ] - probs[2, ]
}
if (design$kMax > 1) {
if (design$sided == 2) {
designCharacteristics$futilityProbabilities <- rep(0, design$kMax - 1)
} else {
designCharacteristics$futilityProbabilities <- probs[1, 1:(design$kMax - 1)]
designCharacteristics$futilityProbabilities[design$futilityBounds == C_FUTILITY_BOUNDS_DEFAULT] <- 0
}
designCharacteristics$power <- .getNoEarlyEfficacyZeroCorrectedValues(
design, designCharacteristics$power
)
designCharacteristics$rejectionProbabilities <- .getNoEarlyEfficacyZeroCorrectedValues(
design, designCharacteristics$rejectionProbabilities
)
}
designCharacteristics$information <- informationRates * shift
designCharacteristics$averageSampleNumber1 <- .getAverageSampleNumber(
design$kMax,
design$informationRates, probs, shift, nFixed
)
if (design$sided == 2) {
decisionMatrix <- matrix(c(-design$criticalValues, design$criticalValues), nrow = 2, byrow = TRUE)
} else {
decisionMatrix <- matrix(c(
design$futilityBounds, C_FUTILITY_BOUNDS_DEFAULT,
design$criticalValues
), nrow = 2, byrow = TRUE)
}
probs0 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber0 <- .getAverageSampleNumber(
design$kMax,
design$informationRates, probs0, shift, nFixed
)
if (design$sided == 2) {
decisionMatrix <- matrix(c(
-design$criticalValues - sqrt(shift * informationRates) / 2,
design$criticalValues - sqrt(shift * informationRates) / 2
), nrow = 2, byrow = TRUE)
} else {
shiftedFutilityBounds <- design$futilityBounds - sqrt(shift * informationRates[1:(design$kMax - 1)]) / 2
shiftedFutilityBounds[design$futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <- C_FUTILITY_BOUNDS_DEFAULT
decisionMatrix <- matrix(c(shiftedFutilityBounds, C_FUTILITY_BOUNDS_DEFAULT, design$criticalValues -
sqrt(shift * informationRates) / 2), nrow = 2, byrow = TRUE)
}
probs01 <- .getGroupSequentialProbabilities(decisionMatrix, informationRates)
designCharacteristics$averageSampleNumber01 <- .getAverageSampleNumber(
design$kMax, design$informationRates, probs01, shift, nFixed
)
}
design$criticalValues[design$criticalValues >= 7.5 - 1e-8] <- Inf
designCharacteristics$.setParameterType("shift", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("power", C_PARAM_GENERATED)
designCharacteristics$.setParameterType(".probs", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("rejectionProbabilities", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("information", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("futilityProbabilities", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("averageSampleNumber0", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("averageSampleNumber01", C_PARAM_GENERATED)
designCharacteristics$.setParameterType("averageSampleNumber1", C_PARAM_GENERATED)
designCharacteristics$inflationFactor <- shift / nFixed
designCharacteristics$.setParameterType("inflationFactor", C_PARAM_GENERATED)
if (is.na(designCharacteristics$inflationFactor) ||
designCharacteristics$inflationFactor > 4 || designCharacteristics$inflationFactor < 1 - 1e-05) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "inflation factor cannot be calculated")
}
return(designCharacteristics)
}
.getAverageSampleNumber <- function(kMax, informationRates, probs, shift, nFixed) {
if (nrow(probs) == 3) {
return((shift - sum((probs[3, 1:(kMax - 1)] - probs[2, 1:(kMax - 1)] + probs[1, 1:(kMax - 1)]) *
(informationRates[kMax] - informationRates[1:(kMax - 1)]) * shift)) / nFixed)
} else {
return((shift - sum((probs[5, 1:(kMax - 1)] -
probs[4, 1:(kMax - 1)] + probs[3, 1:(kMax - 1)] - probs[2, 1:(kMax - 1)] + probs[1, 1:(kMax - 1)]) *
(informationRates[kMax] - informationRates[1:(kMax - 1)]) * shift)) / nFixed)
}
}
#'
#' @title
#' Get Power And Average Sample Number
#'
#' @description
#' Returns the power and average sample number of the specified design.
#'
#' @inheritParams param_design
#' @inheritParams param_theta
#' @inheritParams param_nMax
#'
#' @details
#' This function returns the power and average sample number (ASN) of the specified
#' design for the prototype case which is testing H0: mu = mu0 in a one-sample design.
#' \code{theta} represents the standardized effect \code{(mu - mu0) / sigma} and power and ASN
#' is calculated for maximum sample size \code{nMax}.
#' For other designs than the one-sample test of a mean the standardized effect needs to be adjusted accordingly.
#'
#' @return Returns a \code{\link{PowerAndAverageSampleNumberResult}} object.
#' The following generics (R generic functions) are available for this result object:
#' \itemize{
#' \item \code{\link[=names.FieldSet]{names()}} to obtain the field names,
#' \item \code{\link[=print.FieldSet]{print()}} to print the object,
#' \item \code{\link[=summary.ParameterSet]{summary()}} to display a summary of the object,
#' \item \code{\link[=plot.ParameterSet]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.PowerAndAverageSampleNumberResult]{as.data.frame()}}
#' to coerce the object to a \code{\link[base]{data.frame}},
#' \item \code{\link[=as.matrix.FieldSet]{as.matrix()}} to coerce the object to a \code{\link[base]{matrix}}.
#' }
#' @template how_to_get_help_for_generics
#'
#' @family design functions
#'
#' @template examples_get_power_and_average_sample_number
#'
#' @export
#'
getPowerAndAverageSampleNumber <- function(design, theta = seq(-1, 1, 0.02), nMax = 100) {
.assertIsTrialDesign(design)
.assertIsSingleNumber(nMax, "nMax")
.assertIsInClosedInterval(nMax, "nMax", lower = 1, upper = NULL)
return(PowerAndAverageSampleNumberResult(design = design, theta = theta, nMax = nMax))
}
.getSimulatedRejectionsDelayedResponse <- function(delta, informationRates, delayedInformation,
contRegionUpper, contRegionLower, decisionCriticalValues, iterations, seed = NA_real_) {
seed <- .setSeed(seed)
kMax <- length(informationRates)
zVector <- numeric(kMax)
reject <- 0L
for (i in 1:iterations) {
for (stage in 1:kMax) {
if (stage == 1) {
zVector[stage] <- stats::rnorm(1, delta * sqrt(informationRates[1]), 1)
} else {
zVector[stage] <- (sqrt(informationRates[stage - 1]) * zVector[stage - 1] +
sqrt(informationRates[stage] - informationRates[stage - 1]) *
stats::rnorm(1, delta * sqrt(informationRates[stage] - informationRates[stage - 1]), 1)) /
sqrt(informationRates[stage])
}
if (!is.na(decisionCriticalValues[stage]) && stage < kMax &&
(zVector[stage] > contRegionUpper[stage] || zVector[stage] < contRegionLower[stage])) {
if ((sqrt(informationRates[stage]) * zVector[stage] + sqrt(delayedInformation[stage]) *
stats::rnorm(1, delta * sqrt(delayedInformation[stage]), 1)) /
sqrt(informationRates[stage] + delayedInformation[stage]) > decisionCriticalValues[stage]) {
reject <- reject + 1L
}
break
}
if (stage == kMax && zVector[stage] > decisionCriticalValues[stage]) {
reject <- reject + 1L
}
}
}
simulatedAlpha <- reject / iterations
return(list(
simulatedAlpha = simulatedAlpha,
delta = delta,
iterations = iterations,
seed = seed
))
}
#'
#' @title
#' Simulated Rejections Delayed Response
#'
#' @description
#' Simulates the rejection probability of a delayed response group sequential design with specified parameters.
#'
#' @inheritParams param_design
#' @param delta The delay value.
#' @param iterations The number of simulation iterations.
#' @inheritParams param_seed
#'
#' @details
#' By default, delta = 0, i.e., the Type error rate is simulated.
#'
#' @return Returns a list summarizing the rejection probabilities.
#'
#' @noRd
#'
getSimulatedRejectionsDelayedResponse <- function(design, ..., delta = 0, iterations = 10000, seed = NA_real_) {
.assertIsTrialDesignInverseNormalOrGroupSequential(design)
.assertIsSingleNumber(delta, "delta")
.assertIsValidIterationsAndSeed(iterations, seed, zeroIterationsAllowed = FALSE)
if (!design$.isDelayedResponseDesign()) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'design' must be a delayed response design with specified 'delayedInformation'")
}
startTime <- Sys.time()
result <- .getSimulatedRejectionsDelayedResponse(
delta = delta,
informationRates = design$informationRates,
delayedInformation = design$delayedInformation,
contRegionUpper = design$criticalValues,
contRegionLower = design$futilityBounds,
decisionCriticalValues = design$decisionCriticalValues,
iterations = iterations,
seed = seed
)
simulatedAlpha <- result$simulatedAlpha
stdError <- sqrt(simulatedAlpha * (1 - simulatedAlpha) / iterations)
ciLower <- simulatedAlpha - 1.96 * stdError
ciUpper <- simulatedAlpha + 1.96 * stdError
result$confidenceIntervall <- c(ciLower, ciUpper)
# simulated Type I error rate is within the 95% error bounds
result$alphaWithin95ConfidenceIntervall <- design$alpha > ciLower && design$alpha < ciUpper
result$time <- Sys.time() - startTime
return(result)
}
Try the rpact package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.