Nothing
R/f_design_plan_survival.R
In rpact: Confirmatory Adaptive Clinical Trial Design and Analysis
Defines functions
getPowerSurvival
getSampleSizeSurvival
getNumberOfSubjects
getEventProbabilities
.getNumberOfSubjects
.getSampleSizeSequentialSurvival
.getSampleSizeFixedSurvival
.initDesignPlanSurvival
.initDesignPlanSurvivalByPiecewiseSurvivalTimeObject
.createDesignPlanSurvival
.getSampleSizeSurvival
.calculateSampleSizeSurvival
.getEffectScaleBoundaryDataSurvival
.getEventProbabilities
.getEventsFixed
.isUserDefinedMaxNumberOfSubjects
.getEventProbabilitiesOverall
.getEventProbabilitiesGroupwise
.getEventProbabilityFunctionVec
.getEventProbabilityFunction
.getLambda
.addStudyDurationToDesignPlan
Documented in
getEventProbabilities
getNumberOfSubjects
getPowerSurvival
getSampleSizeSurvival
## |
## | *Sample size and power survival*
## |
## | 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: 7662 $
## | Last changed: $Date: 2024-02-23 12:42:26 +0100 (Fr, 23 Feb 2024) $
## | Last changed by: $Author: pahlke $
## |
#' @include f_core_utilities.R
NULL
.addStudyDurationToDesignPlan <- function(designPlan) {
if (!designPlan$accountForObservationTimes) {
return(invisible())
}
kMax <- designPlan$.design$kMax
if (kMax == 1) {
designPlan$studyDuration <- designPlan$analysisTime[1, ]
designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
designPlan$maxStudyDuration <- designPlan$studyDuration
} else {
designPlan$maxStudyDuration <- designPlan$analysisTime[kMax, ]
designPlan$.setParameterType("maxStudyDuration", C_PARAM_GENERATED)
}
}
.getLambda <- function(..., groupNumber, lambda2, lambda1, hazardRatio, kappa) {
if (groupNumber == 1) {
if (!any(is.na(lambda1))) {
return(lambda1)
}
lambda2 <- lambda2 * hazardRatio^(1 / kappa)
}
return(lambda2)
}
.getEventProbabilityFunction <- function(..., time, piecewiseLambda, piecewiseSurvivalTime, phi, kappa) {
if (length(piecewiseLambda) == 1) {
if (kappa != 1 && phi > 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "Weibull distribution cannot ",
"be used together with specified dropout rate (use simulation instead)",
call. = FALSE
)
}
return(piecewiseLambda / (piecewiseLambda + phi) *
pweibull(time, shape = kappa, scale = 1 / (piecewiseLambda + phi), lower.tail = TRUE, log.p = FALSE))
}
if (length(piecewiseSurvivalTime) != length(piecewiseLambda)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of 'piecewiseSurvivalTime' (", .arrayToString(piecewiseSurvivalTime),
") must be equal to length of 'piecewiseLambda' (", .arrayToString(piecewiseLambda), ")"
)
}
piecewiseSurvivalTime <- .getPiecewiseExpStartTimesWithoutLeadingZero(piecewiseSurvivalTime)
if (kappa != 1) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"Weibull distribution cannot be used for piecewise survival definition",
call. = FALSE
)
}
len <- length(piecewiseSurvivalTime)
for (i in 1:len) {
if (i == 1) {
if (time <= piecewiseSurvivalTime[1]) {
return(piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
(1 - exp(-((piecewiseLambda[1] + phi) * time))))
}
} else if (i == 2) {
cdfPart <- piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
(1 - exp(-((piecewiseLambda[1] + phi) * piecewiseSurvivalTime[1])))
if (time <= piecewiseSurvivalTime[2]) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
cdf <- cdfPart + piecewiseLambda[2] / (piecewiseLambda[2] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[1]) - exp(-piecewiseLambda[2] *
(time - piecewiseSurvivalTime[1]) - phi * time))
return(cdf)
}
} else if (i == 3) {
cdfPart <- cdfPart + piecewiseLambda[2] / (piecewiseLambda[2] + phi) *
exp(-piecewiseLambda[1] * piecewiseSurvivalTime[1]) * (
exp(-phi * piecewiseSurvivalTime[1]) - exp(-piecewiseLambda[2] *
(piecewiseSurvivalTime[2] - piecewiseSurvivalTime[1]) - phi * piecewiseSurvivalTime[2]))
if (time <= piecewiseSurvivalTime[3]) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
piecewiseLambda[2] * (piecewiseSurvivalTime[2] - piecewiseSurvivalTime[1])
cdf <- cdfPart + piecewiseLambda[3] / (piecewiseLambda[3] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[2]) - exp(-piecewiseLambda[3] *
(time - piecewiseSurvivalTime[2]) - phi * time))
return(cdf)
}
} else if (i > 3) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
sum(piecewiseLambda[2:(i - 2)] * (piecewiseSurvivalTime[2:(i - 2)] -
piecewiseSurvivalTime[1:(i - 3)]))
cdfPart <- cdfPart + piecewiseLambda[i - 1] / (piecewiseLambda[i - 1] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[i - 2]) - exp(-piecewiseLambda[i - 1] *
(piecewiseSurvivalTime[i - 1] - piecewiseSurvivalTime[i - 2]) - phi * piecewiseSurvivalTime[i - 1]))
if (time <= piecewiseSurvivalTime[i]) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
sum(piecewiseLambda[2:(i - 1)] * (piecewiseSurvivalTime[2:(i - 1)] - piecewiseSurvivalTime[1:(i - 2)]))
cdf <- cdfPart + piecewiseLambda[i] / (piecewiseLambda[i] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[i - 1]) - exp(-piecewiseLambda[i] *
(time - piecewiseSurvivalTime[i - 1]) - phi * time))
return(cdf)
}
}
}
if (len == 1) {
cdfPart <- piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
(1 - exp(-((piecewiseLambda[1] + phi) * piecewiseSurvivalTime[1])))
} else if (len == 2) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
cdfPart <- cdfPart + piecewiseLambda[len] / (piecewiseLambda[len] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[len - 1]) - exp(-piecewiseLambda[len] *
(piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1]) - phi * piecewiseSurvivalTime[len]))
} else {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
sum(piecewiseLambda[2:(len - 1)] * (piecewiseSurvivalTime[2:(len - 1)] - piecewiseSurvivalTime[1:(len - 2)]))
cdfPart <- cdfPart + piecewiseLambda[len] / (piecewiseLambda[len] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[len - 1]) - exp(-piecewiseLambda[len] *
(piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1]) - phi * piecewiseSurvivalTime[len]))
}
if (len == 1) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
} else {
cdfFactor <- cdfFactor + piecewiseLambda[len] * (piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1])
}
cdf <- cdfPart + piecewiseLambda[len + 1] / (piecewiseLambda[len + 1] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[len]) - exp(-piecewiseLambda[len + 1] *
(time - piecewiseSurvivalTime[len]) - phi * time))
return(cdf)
}
.getEventProbabilityFunctionVec <- function(..., timeVector, piecewiseLambda, piecewiseSurvivalTime, phi, kappa) {
result <- c()
for (time in timeVector) {
result <- c(result, .getEventProbabilityFunction(
time = time, piecewiseLambda = piecewiseLambda,
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa
))
}
return(result)
}
.getEventProbabilitiesGroupwise <- function(
..., time, accrualTimeVector, accrualIntensity, lambda2,
lambda1, piecewiseSurvivalTime, phi, kappa, allocationRatioPlanned, hazardRatio) {
.assertIsSingleNumber(time, "time")
if (length(accrualTimeVector) > 1 && accrualTimeVector[1] == 0) {
accrualTimeVector <- accrualTimeVector[2:length(accrualTimeVector)]
}
accrualTimeVectorLength <- length(accrualTimeVector)
densityIntervals <- accrualTimeVector
if (accrualTimeVectorLength > 1) {
densityIntervals[2:accrualTimeVectorLength] <-
accrualTimeVector[2:accrualTimeVectorLength] -
accrualTimeVector[1:(accrualTimeVectorLength - 1)]
}
if (length(densityIntervals) > 1 && length(accrualIntensity) > 1 &&
length(densityIntervals) != length(accrualIntensity)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'densityIntervals' (", .arrayToString(densityIntervals),
") and 'accrualIntensity' (", .arrayToString(accrualIntensity), ") must have same length"
)
}
densityVector <- accrualIntensity / sum(densityIntervals * accrualIntensity)
eventProbs <- rep(NA_real_, 2)
for (k in 1:accrualTimeVectorLength) {
if (time <= accrualTimeVector[k]) {
for (groupNumber in c(1, 2)) { # two groups: 1 = treatment, 2 = control
lambdaTemp <- .getLambda(
groupNumber = groupNumber,
lambda2 = lambda2,
lambda1 = lambda1,
hazardRatio = hazardRatio,
kappa = kappa
)
inner <- function(x) {
.getEventProbabilityFunctionVec(
timeVector = x, piecewiseLambda = lambdaTemp,
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi[groupNumber], kappa = kappa
)
}
timeValue1 <- 0
if (k > 1) {
timeValue1 <- time - accrualTimeVector[1]
}
eventProbs[groupNumber] <- densityVector[1] * stats::integrate(inner, timeValue1, time)$value
if (k > 2) {
for (j in 2:(k - 1)) {
eventProbs[groupNumber] <- eventProbs[groupNumber] +
densityVector[j] * stats::integrate(
inner, time - accrualTimeVector[j],
time - accrualTimeVector[j - 1]
)$value
}
}
if (k > 1) {
eventProbs[groupNumber] <- eventProbs[groupNumber] +
densityVector[k] * stats::integrate(inner, 0, time - accrualTimeVector[k - 1])$value
}
}
return(eventProbs)
}
}
for (groupNumber in c(1, 2)) {
lambdaTemp <- .getLambda(
groupNumber = groupNumber,
lambda2 = lambda2,
lambda1 = lambda1,
hazardRatio = hazardRatio,
kappa = kappa
)
inner <- function(x) {
.getEventProbabilityFunctionVec(
timeVector = x, piecewiseLambda = lambdaTemp,
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi[groupNumber], kappa = kappa
)
}
eventProbs[groupNumber] <- densityVector[1] *
stats::integrate(inner, time - accrualTimeVector[1], time)$value
if (accrualTimeVectorLength > 1) {
for (j in (2:accrualTimeVectorLength)) {
eventProbs[groupNumber] <- eventProbs[groupNumber] +
densityVector[j] * stats::integrate(
inner, time - accrualTimeVector[j],
time - accrualTimeVector[j - 1]
)$value
}
}
}
return(eventProbs)
}
.getEventProbabilitiesOverall <- function(eventProbs, allocationRatioPlanned) {
return((allocationRatioPlanned * eventProbs[1] + eventProbs[2]) / (1 + allocationRatioPlanned))
}
.isUserDefinedMaxNumberOfSubjects <- function(designPlan) {
if (!is.null(designPlan) && length(designPlan$.getParameterType("maxNumberOfSubjects")) > 0) {
if (designPlan$.getParameterType("maxNumberOfSubjects") == C_PARAM_USER_DEFINED) {
return(TRUE)
}
}
return(!is.null(designPlan$maxNumberOfSubjects) && length(designPlan$maxNumberOfSubjects) == 1 &&
!is.na(designPlan$maxNumberOfSubjects) && designPlan$maxNumberOfSubjects > 0)
}
.getEventsFixed <- function(..., typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
twoSidedPower, alpha, beta, sided, hazardRatio, thetaH0, allocationRatioPlanned) {
typeOfComputation <- match.arg(typeOfComputation)
if (typeOfComputation == "Schoenfeld") {
eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 /
(log(hazardRatio) - log(thetaH0))^2 *
(1 + allocationRatioPlanned)^2 / allocationRatioPlanned
if (twoSidedPower && (sided == 2)) {
up <- 2 * eventsFixed
eventsFixed <- .getOneDimensionalRoot(
function(n) {
return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
(log(hazardRatio) - log(thetaH0)) * sqrt(allocationRatioPlanned) /
(1 + allocationRatioPlanned)) -
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
(log(hazardRatio) - log(thetaH0)) * sqrt(allocationRatioPlanned) /
(1 + allocationRatioPlanned)) - beta)
},
lower = 0.001, upper = up, tolerance = 1e-04,
callingFunctionInformation = ".getEventsFixed"
)
}
return(eventsFixed)
}
if (typeOfComputation == "Freedman") {
eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 *
(1 + hazardRatio * allocationRatioPlanned)^2 / (1 - hazardRatio)^2 /
allocationRatioPlanned
if (twoSidedPower && (sided == 2)) {
up <- 2 * eventsFixed
eventsFixed <- .getOneDimensionalRoot(
function(n) {
return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
sqrt(allocationRatioPlanned) * (1 - hazardRatio) /
(1 + allocationRatioPlanned * hazardRatio)) -
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
sqrt(allocationRatioPlanned) * (1 - hazardRatio) /
(1 + allocationRatioPlanned * hazardRatio)) - beta)
},
lower = 0.001, upper = up, tolerance = 1e-04,
callingFunctionInformation = ".getEventsFixed"
)
}
return(eventsFixed)
}
if (typeOfComputation == "HsiehFreedman") {
eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 *
(1 + hazardRatio)^2 / (1 - hazardRatio)^2 *
(1 + allocationRatioPlanned)^2 / (4 * allocationRatioPlanned)
if (twoSidedPower && sided == 2) {
up <- 2 * eventsFixed
eventsFixed <- .getOneDimensionalRoot(
function(n) {
return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
2 * sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(1 - hazardRatio) / (1 + hazardRatio)) -
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
2 * sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(1 - hazardRatio) / (1 + hazardRatio)) - beta)
},
lower = 0.001, upper = up, tolerance = 1e-04,
callingFunctionInformation = ".getEventsFixed"
)
}
return(eventsFixed)
}
}
.getEventProbabilities <- function(..., time, accrualTimeVector, accrualIntensity, lambda2,
lambda1, piecewiseSurvivalTime, phi, kappa, allocationRatioPlanned, hazardRatio) {
eventProbs <- .getEventProbabilitiesGroupwise(
time = time, accrualTimeVector = accrualTimeVector,
accrualIntensity = accrualIntensity, lambda2 = lambda2,
lambda1 = lambda1, piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio
)
return(.getEventProbabilitiesOverall(eventProbs, allocationRatioPlanned))
}
.getEffectScaleBoundaryDataSurvival <- function(designPlan) {
design <- designPlan$.design
thetaH0 <- designPlan$thetaH0
cumulativeEventsPerStage <- designPlan$cumulativeEventsPerStage
allocationRatioPlanned <- designPlan$allocationRatioPlanned
directionUpper <- designPlan$directionUpper
if (design$kMax == 1) {
nParameters <- length(cumulativeEventsPerStage)
} else {
nParameters <- ncol(cumulativeEventsPerStage)
}
directionUpper[is.na(directionUpper)] <- TRUE
if (length(allocationRatioPlanned) == 1) {
allocationRatioPlanned <- rep(allocationRatioPlanned, nParameters)
}
futilityBounds <- design$futilityBounds
futilityBounds[!is.na(futilityBounds) & futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <- NA_real_
criticalValues <- design$criticalValues
criticalValuesEffectScaleUpper <- matrix(, nrow = design$kMax, ncol = nParameters)
criticalValuesEffectScaleLower <- matrix(, nrow = design$kMax, ncol = nParameters)
futilityBoundsEffectScaleUpper <- matrix(, nrow = design$kMax - 1, ncol = nParameters)
futilityBoundsEffectScaleLower <- matrix(, nrow = design$kMax - 1, ncol = nParameters)
for (j in (1:nParameters)) {
if (design$sided == 1) {
criticalValuesEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * criticalValues *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[, j])))
} else {
criticalValuesEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * criticalValues *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[, j])))
criticalValuesEffectScaleLower[, j] <- thetaH0 * (exp(-(2 * directionUpper[j] - 1) * criticalValues *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[, j])))
}
if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
futilityBoundsEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * futilityBounds *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[1:(design$kMax - 1), j])))
}
if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
(design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
futilityBoundsEffectScaleLower[, j] <- thetaH0 * (exp(-(2 * directionUpper[j] - 1) * futilityBounds *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[1:(design$kMax - 1), j])))
}
}
return(list(
criticalValuesEffectScaleUpper = matrix(criticalValuesEffectScaleUpper, nrow = design$kMax),
criticalValuesEffectScaleLower = matrix(criticalValuesEffectScaleLower, nrow = design$kMax),
futilityBoundsEffectScaleUpper = matrix(futilityBoundsEffectScaleUpper, nrow = design$kMax - 1),
futilityBoundsEffectScaleLower = matrix(futilityBoundsEffectScaleLower, nrow = design$kMax - 1)
))
}
.calculateSampleSizeSurvival <- function(designPlan) {
# Fixed
designPlan <- .getSampleSizeFixedSurvival(designPlan)
# Sequential
if (designPlan$.design$kMax > 1) {
designCharacteristics <- getDesignCharacteristics(designPlan$.design)
designPlan <- .getSampleSizeSequentialSurvival(designPlan, designCharacteristics)
}
if (designPlan$accountForObservationTimes && !any(is.na(designPlan$followUpTime)) &&
all(designPlan$followUpTime == C_FOLLOW_UP_TIME_DEFAULT)) {
designPlan$.setParameterType("followUpTime", C_PARAM_DEFAULT_VALUE)
}
.addEffectScaleBoundaryDataToDesignPlan(designPlan)
designPlan$singleEventsPerStage <- .getSingleEventsPerStage(designPlan$cumulativeEventsPerStage)
.addDeprecatedFieldValues(designPlan, "eventsPerStage", designPlan$cumulativeEventsPerStage)
paramName <- ifelse(designPlan$.design$kMax == 1, "nFixed", "maxNumberOfSubjects")
if (designPlan$.getParameterType(paramName) == C_PARAM_GENERATED &&
designPlan$.accrualTime$.getParameterType("maxNumberOfSubjects") != C_PARAM_GENERATED &&
all(designPlan$accrualIntensity < 1)) {
numberOfDefinedAccrualIntensities <- length(designPlan$accrualIntensity)
accrualTime <- designPlan$accrualTime
if (length(accrualTime) > 0 && accrualTime[1] != 0) {
accrualTime <- c(0, accrualTime)
}
if (any(designPlan$accrualIntensity < 0)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'accrualIntensityRelative' (",
.arrayToString(designPlan$accrualIntensity), ") must be >= 0"
)
}
designPlan$accrualIntensityRelative <- designPlan$accrualIntensity
if (identical(designPlan$accrualIntensityRelative, C_ACCRUAL_INTENSITY_DEFAULT)) {
designPlan$.setParameterType("accrualIntensityRelative", C_PARAM_NOT_APPLICABLE)
} else {
designPlan$.setParameterType(
"accrualIntensityRelative",
designPlan$.getParameterType("accrualIntensity")
)
}
accrualIntensityAbsolute <- numeric()
for (maxNumberOfSubjects in designPlan[[paramName]]) {
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = designPlan$accrualIntensityRelative,
accrualIntensityType = "relative",
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualIntensityAbsolute <- c(accrualIntensityAbsolute, accrualSetup$accrualIntensity)
}
designPlan$accrualIntensity <- accrualIntensityAbsolute
designPlan$.setParameterType("accrualIntensity", C_PARAM_GENERATED)
if (numberOfDefinedAccrualIntensities > 1) {
paramName <- NULL
if (designPlan$.getParameterType("pi1") == C_PARAM_USER_DEFINED ||
designPlan$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE ||
designPlan$.getParameterType("pi2") == C_PARAM_USER_DEFINED) {
paramName <- "pi1"
} else if (designPlan$.getParameterType("median1") == C_PARAM_USER_DEFINED ||
designPlan$.getParameterType("median2") == C_PARAM_USER_DEFINED) {
paramName <- "median1"
}
if (!is.null(paramName)) {
paramValue <- designPlan[[paramName]]
if (!is.null(paramValue) && length(paramValue) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the definition of relative accrual intensities ",
"(all 'accrualIntensity' values < 1) ",
"is only available for a single value ",
"(", paramName, " = ", .arrayToString(
paramValue,
vectorLookAndFeelEnabled = TRUE
), ")"
)
}
}
}
}
designPlan$maxNumberOfEvents <- designPlan$cumulativeEventsPerStage[designPlan$.design$kMax, ]
designPlan$.setParameterType("maxNumberOfEvents", C_PARAM_GENERATED)
if (!any(is.na(designPlan$followUpTime))) {
if (any(designPlan$followUpTime < -1e-02)) {
warning("Accrual duration longer than maximal study ",
"duration (time to maximal number of events); followUpTime = ",
.arrayToString(designPlan$followUpTime),
call. = FALSE
)
}
} else {
warning("Follow-up time could not be calculated for hazardRatio = ",
.arrayToString(designPlan$hazardRatio[indices]),
call. = FALSE
)
}
if (designPlan$.getParameterType("accountForObservationTimes") != C_PARAM_USER_DEFINED) {
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_NOT_APPLICABLE)
}
designPlan$.setParameterType("chi", C_PARAM_NOT_APPLICABLE)
.addStudyDurationToDesignPlan(designPlan)
return(designPlan)
}
# Hidden parameter:
# @param accountForObservationTimes If \code{accountForObservationTimes = TRUE}, the number of
# subjects is calculated assuming specific accrual and follow-up time, default is \code{TRUE}
# (see details).
# If \code{accountForObservationTimes = FALSE}, only the event rates are used for the calculation
# of the maximum number of subjects.
# \code{accountForObservationTimes} can be selected as \code{FALSE}. In this case,
# the number of subjects is calculated from the event probabilities only.
# This kind of computation does not account for the specific accrual pattern and survival distribution.
.getSampleSizeSurvival <- function(...,
design = NULL,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0 = 1,
pi2 = NA_real_,
pi1 = NA_real_,
allocationRatioPlanned = NA_real_,
accountForObservationTimes = TRUE,
eventTime = C_EVENT_TIME_DEFAULT,
accrualTime = C_ACCRUAL_TIME_DEFAULT,
accrualIntensity = C_ACCRUAL_INTENSITY_DEFAULT,
accrualIntensityType = c("auto", "absolute", "relative"),
kappa = 1,
piecewiseSurvivalTime = NA_real_,
lambda2 = NA_real_,
lambda1 = NA_real_,
median1 = NA_real_,
median2 = NA_real_,
followUpTime = NA_real_,
maxNumberOfSubjects = NA_real_,
dropoutRate1 = 0,
dropoutRate2 = dropoutRate1,
dropoutTime = NA_real_,
hazardRatio = NA_real_,
userFunctionCallEnabled = TRUE) {
designPlan <- .createDesignPlanSurvival(
objectType = "sampleSize",
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0,
pi2 = pi2,
pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2,
lambda1 = lambda1,
median1 = median1,
median2 = median2,
followUpTime = followUpTime,
maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = userFunctionCallEnabled
)
return(.calculateSampleSizeSurvival(designPlan))
}
.createDesignPlanSurvival <- function(..., objectType = c("sampleSize", "power"),
design,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0, pi2, pi1,
allocationRatioPlanned,
accountForObservationTimes,
eventTime,
accrualTime,
accrualIntensity,
accrualIntensityType,
kappa,
piecewiseSurvivalTime,
lambda2,
lambda1,
median1,
median2,
followUpTime = NA_real_,
directionUpper = NA,
maxNumberOfEvents = NA_real_,
maxNumberOfSubjects,
dropoutRate1,
dropoutRate2,
dropoutTime,
hazardRatio,
userFunctionCallEnabled = TRUE) {
objectType <- match.arg(objectType)
typeOfComputation <- .matchArgument(typeOfComputation, "Schoenfeld")
if (userFunctionCallEnabled) {
.assertIsTrialDesignInverseNormalOrGroupSequential(design)
.assertIsValidAlphaAndBeta(design$alpha, design$beta)
.assertIsValidSidedParameter(design$sided)
.assertIsSingleLogical(accountForObservationTimes, "accountForObservationTimes", naAllowed = TRUE)
.assertIsSingleNumber(thetaH0, "thetaH0")
.assertIsValidThetaH0(thetaH0, endpoint = "survival", groups = 2)
.assertIsValidKappa(kappa)
directionUpper <- .assertIsValidDirectionUpper(directionUpper,
design$sided, objectType,
userFunctionCallEnabled = TRUE
)
if (objectType == "power") {
.assertIsSingleNumber(maxNumberOfEvents, "maxNumberOfEvents")
.assertIsInClosedInterval(maxNumberOfEvents, "maxNumberOfEvents",
lower = 1, upper = maxNumberOfSubjects
)
}
if (!any(is.na(pi1)) && (any(pi1 <= 0) || any(pi1 >= 1))) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"event rate 'pi1' (", .arrayToString(pi1), ") is out of bounds (0; 1)"
)
}
if (!any(is.na(pi2)) && (any(pi2 <= 0) || any(pi2 >= 1))) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"event rate 'pi2' (", .arrayToString(pi2), ") is out of bounds (0; 1)"
)
}
if (design$sided == 2 && thetaH0 != 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"two-sided case is implemented for superiority testing only (i.e., thetaH0 = 1)"
)
}
if (thetaH0 <= 0) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"null hypothesis hazard ratio is not allowed be negative or zero"
)
}
if (!(typeOfComputation %in% c("Schoenfeld", "Freedman", "HsiehFreedman"))) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"computation type ('", typeOfComputation, "') must be one of the following: ",
"'Schoenfeld', 'Freedman', or 'HsiehFreedman' "
)
}
if (typeOfComputation != "Schoenfeld" && thetaH0 != 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"Freedman test calculation is possible only for superiority testing (thetaH0 != 1)"
)
}
}
if (is.numeric(accrualTime) && all(is.na(accrualTime))) {
accrualTime <- C_ACCRUAL_TIME_DEFAULT
}
if (all(is.na(accrualIntensity))) {
accrualIntensity <- C_ACCRUAL_INTENSITY_DEFAULT
}
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualSetup$.validate()
if (is.na(allocationRatioPlanned)) {
allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
}
if (userFunctionCallEnabled) {
.assertIsValidAllocationRatioPlannedSampleSize(allocationRatioPlanned, maxNumberOfSubjects)
}
designPlan <- TrialDesignPlanSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime,
accrualTime = accrualSetup$.getAccrualTimeWithoutLeadingZero(),
accrualIntensity = accrualSetup$accrualIntensity,
kappa = kappa,
followUpTime = followUpTime,
maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)
.setValueAndParameterType(
designPlan, "allocationRatioPlanned",
allocationRatioPlanned, C_ALLOCATION_RATIO_DEFAULT
)
.setValueAndParameterType(designPlan, "dropoutRate1", dropoutRate1, C_DROP_OUT_RATE_1_DEFAULT)
.setValueAndParameterType(designPlan, "dropoutRate2", dropoutRate2, C_DROP_OUT_RATE_2_DEFAULT)
.setValueAndParameterType(designPlan, "dropoutTime", dropoutTime, C_DROP_OUT_TIME_DEFAULT)
.setValueAndParameterType(designPlan, "kappa", kappa, 1)
designPlan$.setObjectType(objectType)
designPlan$criticalValuesPValueScale <- matrix(design$stageLevels, ncol = 1)
if (design$sided == 2) {
designPlan$criticalValuesPValueScale <- designPlan$criticalValuesPValueScale * 2
designPlan$.setParameterType("criticalValuesPValueScale", C_PARAM_GENERATED)
}
if (any(design$futilityBounds > C_FUTILITY_BOUNDS_DEFAULT)) {
designPlan$futilityBoundsPValueScale <- matrix(1 - stats::pnorm(design$futilityBounds), ncol = 1)
designPlan$.setParameterType("futilityBoundsPValueScale", C_PARAM_GENERATED)
}
designPlan$.accrualTime <- accrualSetup
designPlan$totalAccrualTime <- accrualSetup$accrualTime[length(accrualSetup$accrualTime)]
if (length(accrualSetup$accrualTime) > 2) {
designPlan$.setParameterType("totalAccrualTime", C_PARAM_GENERATED)
} else {
designPlan$.setParameterType("totalAccrualTime", C_PARAM_NOT_APPLICABLE)
}
if (is.na(maxNumberOfSubjects)) {
if (!is.na(accrualSetup$maxNumberOfSubjects)) {
designPlan$maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects
designPlan$.setParameterType(
"maxNumberOfSubjects",
accrualSetup$.getParameterType("maxNumberOfSubjects")
)
}
} else if (maxNumberOfSubjects == 0) {
designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
} else {
designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_USER_DEFINED)
}
if (identical(as.integer(accrualSetup$accrualTime), C_ACCRUAL_TIME_DEFAULT) ||
identical(
as.integer(c(0L, accrualSetup$.getAccrualTimeWithoutLeadingZero())),
C_ACCRUAL_TIME_DEFAULT
)) {
designPlan$.setParameterType("accrualTime", C_PARAM_DEFAULT_VALUE)
} else {
designPlan$.setParameterType("accrualTime", accrualSetup$.getParameterType("accrualTime"))
}
if (length(designPlan$accrualIntensity) == 1 &&
designPlan$accrualIntensity == C_ACCRUAL_INTENSITY_DEFAULT) {
designPlan$.setParameterType("accrualIntensity", C_PARAM_DEFAULT_VALUE)
} else {
designPlan$.setParameterType(
"accrualIntensity",
accrualSetup$.getParameterType("accrualIntensity")
)
}
.assertIsSingleNumber(designPlan$eventTime, "eventTime")
.assertIsSingleNumber(designPlan$allocationRatioPlanned, "allocationRatioPlanned")
.assertIsSingleNumber(designPlan$kappa, "kappa")
if (objectType == "power") {
.assertIsValidMaxNumberOfSubjects(designPlan$maxNumberOfSubjects)
}
.assertIsSingleNumber(designPlan$dropoutRate1, "dropoutRate1")
.assertIsSingleNumber(designPlan$dropoutRate2, "dropoutRate2")
.assertIsSingleNumber(designPlan$dropoutTime, "dropoutTime")
if (objectType == "power") {
pi1Default <- C_PI_1_DEFAULT
} else {
pi1Default <- C_PI_1_SAMPLE_SIZE_DEFAULT
}
designPlan$.piecewiseSurvivalTime <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime, lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
hazardRatio = hazardRatio, pi1 = pi1, pi2 = pi2, eventTime = eventTime, kappa = kappa,
.pi1Default = pi1Default
)
designPlan$.setParameterType("kappa", designPlan$.piecewiseSurvivalTime$.getParameterType("kappa"))
if (designPlan$.piecewiseSurvivalTime$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE &&
length(designPlan$.piecewiseSurvivalTime$pi1) > 1 &&
length(accrualSetup$accrualIntensity) > 1 && all(accrualSetup$accrualIntensity < 1)) {
designPlan$.piecewiseSurvivalTime$pi1 <- designPlan$.piecewiseSurvivalTime$pi1[1]
warning("Only the first default 'pi1' (", designPlan$.piecewiseSurvivalTime$pi1, ") was used ",
"because the accrual intensities (", .arrayToString(accrualSetup$accrualIntensity), ") ",
"were defined relative (all accrual intensities are < 1)",
call. = FALSE
)
}
.initDesignPlanSurvival(designPlan)
designPlan$.setParameterType("followUpTime", C_PARAM_NOT_APPLICABLE)
if (designPlan$accountForObservationTimes) {
.assertIsSingleNumber(dropoutRate1, "dropoutRate1")
.assertIsSingleNumber(dropoutRate2, "dropoutRate2")
.assertIsSingleNumber(dropoutTime, "dropoutTime")
if (!is.na(dropoutTime) && dropoutTime <= 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'dropoutTime' (", dropoutTime, ") must be > 0")
}
if (dropoutRate1 < 0 || dropoutRate1 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate1' (", dropoutRate1, ") is out of bounds [0; 1)"
)
}
if (dropoutRate2 < 0 || dropoutRate2 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate2' (", dropoutRate2, ") is out of bounds [0; 1)"
)
}
if (!is.na(eventTime) && eventTime <= 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'eventTime' (", eventTime, ") must be > 0")
}
.assertIsValidAccrualTime(accrualSetup$accrualTime)
.assertIsValidFollowUpTime(followUpTime)
.setValueAndParameterType(designPlan, "followUpTime", followUpTime, C_FOLLOW_UP_TIME_DEFAULT)
if (.isUserDefinedMaxNumberOfSubjects(designPlan) && !is.null(followUpTime) &&
length(followUpTime) == 1 && !is.na(followUpTime)) {
warning("Follow-up time will be calculated, value entered (",
followUpTime, ") is not taken into account",
call. = FALSE
)
} else if (is.na(followUpTime)) {
designPlan$followUpTime <- C_FOLLOW_UP_TIME_DEFAULT
designPlan$.setParameterType("followUpTime", C_PARAM_DEFAULT_VALUE)
}
if (objectType == "power") {
designPlan$followUpTime <- NA_real_
designPlan$.setParameterType("followUpTime", C_PARAM_NOT_APPLICABLE)
}
} else {
for (p in c(
"accrualTime", "accrualIntensity",
"eventTime", "dropoutRate1", "dropoutRate2", "dropoutTime", "followUpTime",
"analysisTime", "studyDuration"
)) {
designPlan$.setParameterType(p, C_PARAM_NOT_APPLICABLE)
}
if (designPlan$.getParameterType("accrualTime") == C_PARAM_USER_DEFINED ||
!identical(accrualTime, C_ACCRUAL_TIME_DEFAULT)) {
designPlan$.warnInCaseArgumentExists(accrualSetup$accrualTime, "accrualTime")
}
designPlan$.warnInCaseArgumentExists(dropoutRate1, "dropoutRate1")
designPlan$.warnInCaseArgumentExists(dropoutRate2, "dropoutRate2")
if (!identical(dropoutTime, C_DROP_OUT_TIME_DEFAULT)) {
designPlan$.warnInCaseArgumentExists(dropoutTime, "dropoutTime")
}
designPlan$.warnInCaseArgumentExists(maxNumberOfSubjects, "maxNumberOfSubjects")
if (!identical(followUpTime, C_FOLLOW_UP_TIME_DEFAULT)) {
designPlan$.warnInCaseArgumentExists(followUpTime, "followUpTime")
}
}
.setValueAndParameterType(designPlan, "directionUpper", directionUpper, TRUE)
if (objectType == "power") {
.setValueAndParameterType(designPlan, "maxNumberOfEvents", maxNumberOfEvents, NA_real_)
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_NOT_APPLICABLE)
}
return(designPlan)
}
.initDesignPlanSurvivalByPiecewiseSurvivalTimeObject <- function(designPlan, pwstSetup) {
designPlan$pi1 <- pwstSetup$pi1
designPlan$.setParameterType("pi1", pwstSetup$.getParameterType("pi1"))
designPlan$pi2 <- pwstSetup$pi2
designPlan$.setParameterType("pi2", pwstSetup$.getParameterType("pi2"))
designPlan$hazardRatio <- pwstSetup$hazardRatio
designPlan$.setParameterType("hazardRatio", pwstSetup$.getParameterType("hazardRatio"))
designPlan$lambda1 <- pwstSetup$lambda1
designPlan$.setParameterType("lambda1", pwstSetup$.getParameterType("lambda1"))
designPlan$lambda2 <- pwstSetup$lambda2
designPlan$.setParameterType("lambda2", pwstSetup$.getParameterType("lambda2"))
designPlan$median1 <- pwstSetup$median1
designPlan$.setParameterType("median1", pwstSetup$.getParameterType("median1"))
designPlan$median2 <- pwstSetup$median2
designPlan$.setParameterType("median2", pwstSetup$.getParameterType("median2"))
designPlan$piecewiseSurvivalTime <- pwstSetup$piecewiseSurvivalTime
designPlan$.setParameterType(
"piecewiseSurvivalTime",
pwstSetup$.getParameterType("piecewiseSurvivalTime")
)
designPlan$eventTime <- pwstSetup$eventTime
designPlan$.setParameterType("eventTime", pwstSetup$.getParameterType("eventTime"))
if (pwstSetup$.isLambdaBased()) {
return(length(designPlan$hazardRatio))
}
return(length(designPlan$pi1))
}
.initDesignPlanSurvival <- function(designPlan) {
numberOfResults <- .initDesignPlanSurvivalByPiecewiseSurvivalTimeObject(
designPlan, designPlan$.piecewiseSurvivalTime
)
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
if (length(designPlan$accountForObservationTimes) == 0 ||
is.na(designPlan$accountForObservationTimes) ||
!designPlan$accountForObservationTimes) {
designPlan$accountForObservationTimes <- TRUE
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_DEFAULT_VALUE)
}
if (!designPlan$accountForObservationTimes) {
designPlan$accountForObservationTimes <- TRUE
warning("'accountForObservationTimes' was set to TRUE ",
"because piecewise exponential survival function is enabled",
call. = FALSE
)
}
} else {
if (.isUserDefinedMaxNumberOfSubjects(designPlan)) {
if (length(designPlan$accountForObservationTimes) != 0 &&
!is.na(designPlan$accountForObservationTimes) &&
!designPlan$accountForObservationTimes) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'accountForObservationTimes' must be TRUE because 'maxNumberOfSubjects' is > 0"
)
}
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_GENERATED)
designPlan$accountForObservationTimes <- TRUE
} else {
if (length(designPlan$accountForObservationTimes) == 0 ||
is.na(designPlan$accountForObservationTimes)) {
designPlan$accountForObservationTimes <- FALSE
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_DEFAULT_VALUE)
} else {
designPlan$.setParameterType(
"accountForObservationTimes",
ifelse(designPlan$accountForObservationTimes,
C_PARAM_DEFAULT_VALUE, C_PARAM_USER_DEFINED
)
)
}
}
}
designPlan$.setParameterType("chi", C_PARAM_NOT_APPLICABLE)
if (designPlan$.isSampleSizeObject()) {
designPlan$.setParameterType("directionUpper", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("maxNumberOfEvents", C_PARAM_NOT_APPLICABLE)
}
return(numberOfResults = numberOfResults)
}
.getSampleSizeFixedSurvival <- function(designPlan) {
alpha <- designPlan$getAlpha()
beta <- designPlan$getBeta()
sided <- designPlan$getSided()
twoSidedPower <- designPlan$getTwoSidedPower()
typeOfComputation <- designPlan$typeOfComputation
thetaH0 <- designPlan$thetaH0
pi1 <- designPlan$pi1
pi2 <- designPlan$pi2
allocationRatioPlanned <- designPlan$allocationRatioPlanned
accountForObservationTimes <- designPlan$accountForObservationTimes
accrualTime <- designPlan$accrualTime
kappa <- designPlan$kappa
piecewiseSurvivalTime <- designPlan$piecewiseSurvivalTime
maxNumberOfSubjects <- designPlan$maxNumberOfSubjects
hazardRatio <- designPlan$hazardRatio
.assertIsValidHazardRatio(hazardRatio, thetaH0)
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
numberOfResults <- length(hazardRatio)
} else {
numberOfResults <- length(pi1)
}
designPlan$eventsFixed <- rep(NA_real_, numberOfResults) # number of events
designPlan$nFixed <- rep(NA_real_, numberOfResults) # number of subjects
designPlan$chi <- rep(NA_real_, numberOfResults) # probability of an event
calculateAllocationRatioPlanned <- FALSE
if (allocationRatioPlanned == 0) {
allocationRatioPlannedVec <- rep(NA_real_, numberOfResults)
calculateAllocationRatioPlanned <- TRUE
designPlan$optimumAllocationRatio <- TRUE
designPlan$.setParameterType("optimumAllocationRatio", C_PARAM_USER_DEFINED)
}
userDefinedMaxNumberOfSubjects <- .isUserDefinedMaxNumberOfSubjects(designPlan)
if (userDefinedMaxNumberOfSubjects && allocationRatioPlanned == 0) {
stop(
C_EXCEPTION_TYPE_CONFLICTING_ARGUMENTS,
"determination of optimum allocation ('allocationRatioPlanned' = 0) not possible ",
"for given 'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ")"
)
}
if (userDefinedMaxNumberOfSubjects) {
timeVector <- rep(NA_real_, numberOfResults)
}
designPlan$.calculateFollowUpTime <- FALSE
lambda1 <- designPlan$lambda1
if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
lambda1 <- rep(NA_real_, numberOfResults)
}
for (i in 1:numberOfResults) {
phi <- -c(
log(1 - designPlan$dropoutRate1),
log(1 - designPlan$dropoutRate2)
) / designPlan$dropoutTime
if (!userDefinedMaxNumberOfSubjects) {
if (calculateAllocationRatioPlanned) {
# allocationRatioPlanned = 0 provides optimum sample size
allocationRatioPlanned <- stats::optimize(function(x) {
numberEvents <- .getEventsFixed(
typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
thetaH0 = thetaH0, allocationRatioPlanned = x
)
if (!accountForObservationTimes) {
probEvent <- (x * pi1[i] + pi2) / (1 + x)
} else {
probEvent <- .getEventProbabilities(
time = accrualTime[length(accrualTime)] + designPlan$followUpTime,
accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2, lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime,
phi = phi, kappa = kappa, allocationRatioPlanned = x,
hazardRatio = hazardRatio[i]
)
}
return(numberEvents / probEvent)
}, interval = c(0, 5), tol = 0.0001)$minimum
allocationRatioPlannedVec[i] <- allocationRatioPlanned
}
designPlan$eventsFixed[i] <- .getEventsFixed(
typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
thetaH0 = thetaH0, allocationRatioPlanned = allocationRatioPlanned
)
if (!accountForObservationTimes) {
designPlan$chi[i] <- (allocationRatioPlanned * pi1[i] + pi2) /
(1 + allocationRatioPlanned)
} else {
designPlan$chi[i] <- .getEventProbabilities(
time = accrualTime[length(accrualTime)] + designPlan$followUpTime,
accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio[i]
)
}
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
designPlan$nFixed[i] <- designPlan$eventsFixed[i] / designPlan$chi[i]
} else {
if (length(maxNumberOfSubjects) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of user defined 'maxNumberOfSubjects' (",
.arrayToString(maxNumberOfSubjects), ") must be 1"
)
}
designPlan$.calculateFollowUpTime <- TRUE
designPlan$eventsFixed[i] <- .getEventsFixed(
typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
thetaH0 = thetaH0, allocationRatioPlanned = allocationRatioPlanned
)
designPlan$nFixed[i] <- maxNumberOfSubjects
if (designPlan$eventsFixed[i] > maxNumberOfSubjects) {
if (length(hazardRatio) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
sprintf(
paste0(
"'maxNumberOfSubjects' (%s) is smaller than the number ",
"of events (%.3f) at index %s (hazard ratio = %.3f)"
),
maxNumberOfSubjects, designPlan$eventsFixed[i], i, hazardRatio[i]
)
)
} else {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
sprintf(
paste0(
"'maxNumberOfSubjects' (%s) is smaller than the number ",
"of events (%.3f)"
),
maxNumberOfSubjects, designPlan$eventsFixed[i]
)
)
}
}
up <- 2
iterate <- 1
while (designPlan$eventsFixed[i] / .getEventProbabilities(
time = up, accrualTimeVector = accrualTime, accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2, lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = hazardRatio[i]
) > maxNumberOfSubjects) {
up <- 2 * up
iterate <- iterate + 1
if (iterate > 50) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the number of subjects is too small to reach maximum number of events ",
"(presumably due to drop-out rates), search algorithm failed"
)
}
}
timeVector[i] <- .getOneDimensionalRoot(
function(x) {
designPlan$eventsFixed[i] / .getEventProbabilities(
time = x, accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity, lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = hazardRatio[i]
) - maxNumberOfSubjects
},
lower = 0, upper = up, tolerance = 1e-06, acceptResultsOutOfTolerance = TRUE,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
if (!is.na(timeVector[i])) {
designPlan$chi[i] <- .getEventProbabilities(
time = timeVector[i],
accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity, lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio[i]
)
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
}
}
}
if (calculateAllocationRatioPlanned) {
allocationRatioPlanned <- allocationRatioPlannedVec
designPlan$allocationRatioPlanned <- allocationRatioPlanned
designPlan$.setParameterType("allocationRatioPlanned", C_PARAM_GENERATED)
}
if (userDefinedMaxNumberOfSubjects) {
designPlan$followUpTime <- timeVector - accrualTime[length(accrualTime)]
designPlan$.setParameterType("followUpTime", C_PARAM_GENERATED)
}
designPlan$nFixed2 <- designPlan$nFixed / (1 + allocationRatioPlanned)
designPlan$nFixed1 <- designPlan$nFixed2 * allocationRatioPlanned
designPlan$numberOfSubjects1 <- matrix(designPlan$nFixed1, nrow = 1)
designPlan$numberOfSubjects2 <- matrix(designPlan$nFixed2, nrow = 1)
if (!designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
eventRatio <- allocationRatioPlanned * pi1 / pi2
} else {
eventRatio <- NA_real_
}
# Fixed
designPlan$hazardRatio <- hazardRatio
designPlan$expectedEventsH1 <- designPlan$eventsFixed
designPlan$maxNumberOfSubjects <- designPlan$nFixed
designPlan$numberOfSubjects <- matrix(designPlan$nFixed, nrow = 1)
designPlan$.setParameterType("eventsFixed", C_PARAM_GENERATED)
designPlan$.setParameterType("nFixed1", C_PARAM_GENERATED)
designPlan$.setParameterType("nFixed2", C_PARAM_GENERATED)
designPlan$.setParameterType("nFixed", C_PARAM_GENERATED)
if (designPlan$accountForObservationTimes) {
designPlan$analysisTime <- matrix(accrualTime[length(accrualTime)] +
designPlan$followUpTime, nrow = 1)
designPlan$.setParameterType("analysisTime", C_PARAM_GENERATED)
}
return(designPlan)
}
# note that fixed sample size must be calculated before on 'designPlan'
.getSampleSizeSequentialSurvival <- function(designPlan, designCharacteristics) {
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
numberOfResults <- length(designPlan$hazardRatio)
} else {
numberOfResults <- length(designPlan$pi1)
}
kMax <- designCharacteristics$.design$kMax
designPlan$cumulativeEventsPerStage <- matrix(NA_real_, kMax, numberOfResults)
analysisTime <- matrix(NA_real_, kMax, numberOfResults)
numberOfSubjects <- matrix(NA_real_, kMax, numberOfResults)
designPlan$expectedEventsH0 <- rep(NA_real_, numberOfResults)
designPlan$expectedEventsH01 <- rep(NA_real_, numberOfResults)
designPlan$expectedEventsH1 <- rep(NA_real_, numberOfResults)
expectedNumberOfSubjectsH1 <- rep(NA_real_, numberOfResults)
studyDuration <- rep(NA_real_, numberOfResults)
designPlan$chi <- rep(NA_real_, numberOfResults)
informationRates <- designCharacteristics$information / designCharacteristics$shift
lambda1 <- designPlan$lambda1
if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
lambda1 <- rep(NA_real_, numberOfResults)
}
if (designPlan$accountForObservationTimes && designPlan$.calculateFollowUpTime) {
designPlan$followUpTime <- rep(NA_real_, numberOfResults)
}
for (i in 1:numberOfResults) {
designPlan$cumulativeEventsPerStage[, i] <- designPlan$eventsFixed[i] * informationRates *
designCharacteristics$inflationFactor
if (!designPlan$accountForObservationTimes) {
if (length(designPlan$allocationRatioPlanned) > 1) {
allocationRatioPlanned <- designPlan$allocationRatioPlanned[i]
} else {
allocationRatioPlanned <- designPlan$allocationRatioPlanned
}
designPlan$chi[i] <- (allocationRatioPlanned * designPlan$pi1[i] + designPlan$pi2) /
(1 + allocationRatioPlanned)
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
numberOfSubjects[kMax, i] <- designPlan$cumulativeEventsPerStage[kMax, i] / designPlan$chi[i]
} else {
phi <- -c(log(1 - designPlan$dropoutRate1), log(1 - designPlan$dropoutRate2)) /
designPlan$dropoutTime
if (designPlan$.calculateFollowUpTime) {
if (designPlan$cumulativeEventsPerStage[kMax, i] > designPlan$maxNumberOfSubjects[i]) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
sprintf(
paste0(
"the number of subjects (%s) is smaller than the number ",
"of events (%s) at stage %s"
),
designPlan$maxNumberOfSubjects[i],
designPlan$cumulativeEventsPerStage[kMax, i], i
)
)
}
up <- 2
iterate <- 1
while (designPlan$cumulativeEventsPerStage[kMax, i] / .getEventProbabilities(
time = up,
accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = designPlan$allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
) > designPlan$maxNumberOfSubjects[i]) {
up <- 2 * up
iterate <- iterate + 1
if (iterate > 50) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the number of subjects is too small to reach maximum number of events ",
"(presumably due to drop-out rates)"
)
}
}
totalTime <- .getOneDimensionalRoot(
function(x) {
designPlan$cumulativeEventsPerStage[kMax, i] / designPlan$maxNumberOfSubjects[i] -
.getEventProbabilities(
time = x, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = designPlan$allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = up, tolerance = 1e-06,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
# analysis times
for (j in 1:kMax) {
analysisTime[j, i] <- .getOneDimensionalRoot(
function(x) {
designPlan$cumulativeEventsPerStage[j, i] / designPlan$maxNumberOfSubjects[i] -
.getEventProbabilities(
time = x, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = designPlan$allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = totalTime, tolerance = 1e-06, acceptResultsOutOfTolerance = TRUE,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
}
analysisTime[kMax, i] <- totalTime
designPlan$followUpTime[i] <- totalTime -
designPlan$accrualTime[length(designPlan$accrualTime)]
numberOfSubjects[, i] <- .getNumberOfSubjects(
time = analysisTime[, i],
accrualTime = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = designPlan$maxNumberOfSubjects[i]
)
} else {
if (length(designPlan$allocationRatioPlanned) > 1) {
allocationRatioPlanned <- designPlan$allocationRatioPlanned[i]
} else {
allocationRatioPlanned <- designPlan$allocationRatioPlanned
}
if (is.na(designPlan$followUpTime)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'followUpTime' must be defined because 'designPlan$.calculateFollowUpTime' = FALSE"
)
}
designPlan$chi[i] <- .getEventProbabilities(
time = designPlan$accrualTime[length(designPlan$accrualTime)] + designPlan$followUpTime,
accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
numberOfSubjects[kMax, i] <- designPlan$cumulativeEventsPerStage[kMax, i] / designPlan$chi[i]
# Analysis times
for (j in 1:(kMax - 1)) {
analysisTime[j, i] <- .getOneDimensionalRoot(
function(x) {
designPlan$cumulativeEventsPerStage[j, i] / numberOfSubjects[kMax, i] -
.getEventProbabilities(
time = x, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = designPlan$accrualTime[length(designPlan$accrualTime)] +
designPlan$followUpTime, tolerance = 1e-06,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
}
analysisTime[kMax, i] <- designPlan$accrualTime[length(designPlan$accrualTime)] +
designPlan$followUpTime
numberOfSubjects[, i] <- .getNumberOfSubjects(
time = analysisTime[, i],
accrualTime = designPlan$accrualTime, accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = numberOfSubjects[kMax, i]
)
}
stoppingProbs <- designCharacteristics$rejectionProbabilities +
c(designCharacteristics$futilityProbabilities, 0)
if (all(is.na(designCharacteristics$futilityProbabilities))) {
warning("Expected number of subjects H1 and study duration H1 ",
"cannot be calculated because the futility probabilities ",
"are not applicable for the specified design",
call. = FALSE
)
}
stoppingProbs[kMax] <- 1 - sum(stoppingProbs[1:(kMax - 1)])
studyDuration[i] <- analysisTime[, i] %*% stoppingProbs
expectedNumberOfSubjectsH1[i] <- numberOfSubjects[, i] %*% stoppingProbs
}
designPlan$expectedEventsH0[i] <- designCharacteristics$averageSampleNumber0 *
designPlan$eventsFixed[i]
designPlan$expectedEventsH01[i] <- designCharacteristics$averageSampleNumber01 *
designPlan$eventsFixed[i]
designPlan$expectedEventsH1[i] <- designCharacteristics$averageSampleNumber1 *
designPlan$eventsFixed[i]
designPlan$.setParameterType("expectedEventsH0", C_PARAM_GENERATED)
designPlan$.setParameterType("expectedEventsH01", C_PARAM_GENERATED)
designPlan$.setParameterType("expectedEventsH1", C_PARAM_GENERATED)
designPlan$numberOfSubjects2 <- numberOfSubjects / (1 + designPlan$allocationRatioPlanned)
designPlan$numberOfSubjects1 <- designPlan$numberOfSubjects2 *
designPlan$allocationRatioPlanned
designPlan$.setParameterType("numberOfSubjects1", C_PARAM_GENERATED)
designPlan$.setParameterType("numberOfSubjects2", C_PARAM_GENERATED)
}
if (!is.null(designCharacteristics$rejectionProbabilities)) {
designPlan$rejectPerStage <- matrix(designCharacteristics$rejectionProbabilities,
nrow = designPlan$.design$kMax
)
designPlan$.setParameterType("rejectPerStage", C_PARAM_GENERATED)
designPlan$earlyStop <- sum(designPlan$rejectPerStage[1:(designPlan$.design$kMax - 1), ])
designPlan$.setParameterType("earlyStop", C_PARAM_GENERATED)
}
if (!is.null(designCharacteristics$futilityProbabilities) &&
any(designPlan$.design$futilityBounds != C_FUTILITY_BOUNDS_DEFAULT)) {
designPlan$futilityPerStage <- matrix(designCharacteristics$futilityProbabilities,
nrow = designPlan$.design$kMax - 1
)
designPlan$.setParameterType("futilityPerStage", C_PARAM_GENERATED)
designPlan$futilityStop <- sum(designPlan$futilityPerStage)
designPlan$.setParameterType("futilityStop", C_PARAM_GENERATED)
designPlan$earlyStop <- designPlan$earlyStop + sum(designPlan$futilityPerStage)
}
designPlan$informationRates <- matrix(informationRates, ncol = 1)
if (!is.matrix(numberOfSubjects)) {
designPlan$numberOfSubjects <- matrix(numberOfSubjects[kMax, ], nrow = 1)
} else {
designPlan$numberOfSubjects <- numberOfSubjects
}
designPlan$maxNumberOfSubjects <- designPlan$numberOfSubjects[kMax, ]
if (designPlan$.getParameterType("maxNumberOfSubjects") == C_PARAM_NOT_APPLICABLE ||
length(designPlan$maxNumberOfSubjects) > 1) {
designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
}
designPlan$maxNumberOfSubjects1 <- .getNumberOfSubjects1(
designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
)
designPlan$maxNumberOfSubjects2 <- .getNumberOfSubjects2(
designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
)
designPlan$.setParameterType("maxNumberOfSubjects1", C_PARAM_GENERATED)
designPlan$.setParameterType("maxNumberOfSubjects2", C_PARAM_GENERATED)
if (ncol(designPlan$informationRates) == 1 &&
identical(designPlan$informationRates[, 1], designPlan$.design$informationRates)) {
designPlan$.setParameterType("informationRates", C_PARAM_NOT_APPLICABLE)
} else {
designPlan$.setParameterType("informationRates", C_PARAM_GENERATED)
}
designPlan$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)
designPlan$.setParameterType("cumulativeEventsPerStage", C_PARAM_GENERATED)
if (designPlan$accountForObservationTimes) {
designPlan$analysisTime <- analysisTime
designPlan$expectedNumberOfSubjectsH1 <- expectedNumberOfSubjectsH1
designPlan$studyDuration <- studyDuration
designPlan$studyDurationH1 <- studyDuration # deprecated
designPlan$.setParameterType("analysisTime", C_PARAM_GENERATED)
designPlan$.setParameterType("expectedNumberOfSubjectsH1", C_PARAM_GENERATED)
designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
}
designPlan$.setParameterType("eventsFixed", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("nFixed1", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("nFixed2", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("nFixed", C_PARAM_NOT_APPLICABLE)
if (all(designPlan$allocationRatioPlanned == 1)) {
designPlan$.setParameterType("numberOfSubjects1", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("numberOfSubjects2", C_PARAM_NOT_APPLICABLE)
}
designPlan$.calculateFollowUpTime <- NA
return(designPlan)
}
.getNumberOfSubjects <- function(..., time, accrualTime, accrualIntensity, maxNumberOfSubjects) {
subjectNumbers <- c()
for (timeValue in time) {
if (is.na(timeValue)) {
return(NA_real_)
}
subjectNumbers <- c(
subjectNumbers,
.getNumberOfSubjectsInner(
timeValue = timeValue, accrualTime = accrualTime,
accrualIntensity = accrualIntensity, maxNumberOfSubjects = maxNumberOfSubjects
)
)
}
return(subjectNumbers)
}
#' @title
#' Get Event Probabilities
#'
#' @description
#' Returns the event probabilities for specified parameters at given time vector.
#'
#' @param time A numeric vector with time values.
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#' the end of accrual at specified \code{accrualIntensity} for the specified
#' number of subjects is determined or \code{accrualIntensity} is calculated
#' at fixed end of accrual.
#' @inheritParams param_three_dots
#'
#' @details
#' The function computes the overall event probabilities in a two treatment groups design.
#' For details of the parameters see \code{\link[=getSampleSizeSurvival]{getSampleSizeSurvival()}}.
#'
#' @return Returns a \code{\link{EventProbabilities}} 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.EventProbabilities]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.ParameterSet]{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
#'
#' @template examples_get_event_probabilities
#'
#' @export
#'
getEventProbabilities <- function(time, ...,
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
kappa = 1,
piecewiseSurvivalTime = NA_real_,
lambda2 = NA_real_,
lambda1 = NA_real_,
allocationRatioPlanned = 1,
hazardRatio = NA_real_,
dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
dropoutTime = 12, # C_DROP_OUT_TIME_DEFAULT
maxNumberOfSubjects = NA_real_) {
.warnInCaseOfUnknownArguments(functionName = "getEventProbabilities", ...)
.assertIsNumericVector(time, "time")
.assertIsValidMaxNumberOfSubjects(maxNumberOfSubjects, naAllowed = TRUE)
.assertIsValidAllocationRatioPlannedSampleSize(allocationRatioPlanned, maxNumberOfSubjects)
.assertIsValidKappa(kappa)
.assertIsSingleNumber(hazardRatio, "hazardRatio", naAllowed = TRUE)
if (!is.na(dropoutTime) && dropoutTime <= 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'dropoutTime' (", dropoutTime, ") must be > 0", call. = FALSE)
}
if (dropoutRate1 < 0 || dropoutRate1 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate1' (", dropoutRate1, ") is out of bounds [0; 1)"
)
}
if (dropoutRate2 < 0 || dropoutRate2 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate2' (", dropoutRate2, ") is out of bounds [0; 1)"
)
}
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualTime <- accrualSetup$.getAccrualTimeWithoutLeadingZero()
accrualIntensity <- accrualSetup$accrualIntensity
maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects
setting <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
hazardRatio = hazardRatio, kappa = kappa,
delayedResponseAllowed = TRUE,
.lambdaBased = TRUE
)
if (!setting$delayedResponseEnabled && length(setting$lambda1) > 1 &&
setting$.getParameterType("lambda1") == C_PARAM_USER_DEFINED) {
warning("Only the first 'lambda1' (", lambda1[1], ") was used to calculate event probabilities", call. = FALSE)
setting <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1[1],
hazardRatio = hazardRatio, kappa = kappa,
delayedResponseAllowed = TRUE,
.lambdaBased = TRUE
)
}
piecewiseSurvivalTime <- setting$piecewiseSurvivalTime
lambda2 <- setting$lambda2
lambda1 <- setting$lambda1
hazardRatio <- setting$hazardRatio
phi <- -log(1 - c(dropoutRate1, dropoutRate2)) / dropoutTime
if (length(accrualTime) != length(accrualIntensity)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of 'accrualTime' (", (length(accrualTime) + 1),
") must be equal to length of 'accrualIntensity' (", length(accrualIntensity), ")"
)
}
if (any(accrualIntensity <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualIntensity' must be > 0")
}
if (any(accrualTime <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualTime' must be > 0")
}
if (kappa != 1 && any(phi > 0)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"for Weibull distribution (kappa != 1) drop-out rates (phi) cannot be specified"
)
}
if (any(phi < 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all drop-out rates (phi) must be >= 0")
}
.assertIsNumericVector(lambda2, "lambda2")
if (any(lambda2 <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all rates (lambda2) must be > 0")
}
eventProbabilities <- EventProbabilities(
.piecewiseSurvivalTime = setting,
.accrualTime = accrualSetup,
time = time,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda1 = lambda1,
lambda2 = lambda2,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = hazardRatio,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
maxNumberOfSubjects = maxNumberOfSubjects
)
eventProbabilities$.setParameterType("time", C_PARAM_USER_DEFINED)
eventProbabilities$.setParameterType(
"accrualTime",
accrualSetup$.getParameterType("accrualTime")
)
eventProbabilities$.setParameterType(
"accrualIntensity",
accrualSetup$.getParameterType("accrualIntensity")
)
eventProbabilities$.setParameterType("kappa", setting$.getParameterType("kappa"))
eventProbabilities$.setParameterType(
"piecewiseSurvivalTime",
setting$.getParameterType("piecewiseSurvivalTime")
)
eventProbabilities$.setParameterType("lambda1", setting$.getParameterType("lambda1"))
eventProbabilities$.setParameterType("lambda2", setting$.getParameterType("lambda2"))
.setValueAndParameterType(eventProbabilities, "allocationRatioPlanned", allocationRatioPlanned, 1)
eventProbabilities$.setParameterType("hazardRatio", setting$.getParameterType("hazardRatio"))
.setValueAndParameterType(eventProbabilities, "dropoutRate1", dropoutRate1, C_DROP_OUT_RATE_1_DEFAULT)
.setValueAndParameterType(eventProbabilities, "dropoutRate2", dropoutRate2, C_DROP_OUT_RATE_2_DEFAULT)
.setValueAndParameterType(eventProbabilities, "dropoutTime", dropoutTime, C_DROP_OUT_TIME_DEFAULT)
eventProbabilities$.setParameterType(
"maxNumberOfSubjects",
accrualSetup$.getParameterType("maxNumberOfSubjects")
)
eventProbabilities$cumulativeEventProbabilities <- numeric(0)
eventProbabilities$overallEventProbabilities <- numeric(0) # deprecated
eventProbabilities$eventProbabilities1 <- numeric(0)
eventProbabilities$eventProbabilities2 <- numeric(0)
for (timeValue in time) {
eventProbs <- .getEventProbabilitiesGroupwise(
time = timeValue,
accrualTimeVector = accrualSetup$.getAccrualTimeWithoutLeadingZero(),
accrualIntensity = accrualSetup$accrualIntensity, lambda2 = lambda2,
lambda1 = lambda1, piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio
)
eventProbabilities$cumulativeEventProbabilities <- c(
eventProbabilities$cumulativeEventProbabilities,
.getEventProbabilitiesOverall(eventProbs, allocationRatioPlanned)
)
eventProbabilities$overallEventProbabilities <-
eventProbabilities$cumulativeEventProbabilities # deprecated
eventProbabilities$eventProbabilities1 <- c(
eventProbabilities$eventProbabilities1,
eventProbs[1]
)
eventProbabilities$eventProbabilities2 <- c(
eventProbabilities$eventProbabilities2,
eventProbs[2]
)
}
eventProbabilities$.setParameterType("cumulativeEventProbabilities", C_PARAM_GENERATED)
eventProbabilities$.setParameterType("eventProbabilities1", C_PARAM_GENERATED)
eventProbabilities$.setParameterType("eventProbabilities2", C_PARAM_GENERATED)
return(eventProbabilities)
}
#' @title
#' Get Number Of Subjects
#'
#' @description
#' Returns the number of recruited subjects at given time vector.
#'
#' @param time A numeric vector with time values.
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#' the end of accrual at specified \code{accrualIntensity} for the specified number of
#' subjects is determined or \code{accrualIntensity} is calculated at fixed end of accrual.
#' @inheritParams param_three_dots
#'
#' @details
#' Calculate number of subjects over time range at given accrual time vector
#' and accrual intensity. Intensity can either be defined in absolute or
#' relative terms (for the latter, \code{maxNumberOfSubjects} needs to be defined)\cr
#' The function is used by \code{\link[=getSampleSizeSurvival]{getSampleSizeSurvival()}}.
#'
#' @return Returns a \code{\link{NumberOfSubjects}} 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.NumberOfSubjects]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.ParameterSet]{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
#'
#' @seealso \code{\link{AccrualTime}} for defining the accrual time.
#'
#' @template examples_get_number_of_subjects
#'
#' @export
#'
getNumberOfSubjects <- function(time, ...,
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
maxNumberOfSubjects = NA_real_) {
.warnInCaseOfUnknownArguments(functionName = "getNumberOfSubjects", ...)
.assertIsNumericVector(time, "time")
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualTime <- accrualSetup$.getAccrualTimeWithoutLeadingZero()
accrualIntensity <- accrualSetup$accrualIntensity
maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects
if (length(accrualTime) != length(accrualIntensity)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of 'accrualTime' (", length(accrualTime),
") must be equal to length of 'accrualIntensity' (", length(accrualIntensity), ")"
)
}
if (any(accrualIntensity < 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualIntensity' must be >= 0")
}
if (all(accrualIntensity < 1)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "at least one value of 'accrualIntensity' must be >= 1")
}
if (any(accrualTime <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualTime' must be > 0")
}
numberOfSubjects <- .getNumberOfSubjects(
time = time, accrualTime = accrualTime,
accrualIntensity = accrualIntensity, maxNumberOfSubjects = maxNumberOfSubjects
)
result <- NumberOfSubjects(
.accrualTime = accrualSetup,
time = time,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
maxNumberOfSubjects = maxNumberOfSubjects,
numberOfSubjects = numberOfSubjects
)
result$.setParameterType("time", C_PARAM_USER_DEFINED)
result$.setParameterType("accrualTime", accrualSetup$.getParameterType("accrualTime"))
result$.setParameterType("accrualIntensity", accrualSetup$.getParameterType("accrualIntensity"))
result$.setParameterType("maxNumberOfSubjects", accrualSetup$.getParameterType("maxNumberOfSubjects"))
result$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)
return(result)
}
#' @title
#' Get Sample Size Survival
#'
#' @description
#' Returns the sample size for testing the hazard ratio in a two treatment groups survival design.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_typeOfComputation
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_thetaH0
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_pi1_survival
#' @inheritParams param_pi2_survival
#' @inheritParams param_median1
#' @inheritParams param_median2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_eventTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param followUpTime The assumed (additional) follow-up time for the study, default is \code{6}.
#' The total study duration is \code{accrualTime + followUpTime}.
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#' the follow-up time for the required number of events is determined.
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the number of events and an estimate for the
#' necessary number of subjects for testing the hazard ratio in a survival design.
#' It also calculates the time when the required events are expected under the given
#' assumptions (exponentially, piecewise exponentially, or Weibull distributed survival times
#' and constant or non-constant piecewise accrual).
#' Additionally, an allocation ratio = \code{n1 / n2} can be specified where \code{n1} and \code{n2} are the number
#' of subjects in the two treatment groups.
#'
#' Optional argument \code{accountForObservationTimes}: if \code{accountForObservationTimes = TRUE}, the number of
#' subjects is calculated assuming specific accrual and follow-up time, default is \code{TRUE}.
#'
#' The formula of Kim & Tsiatis (Biometrics, 1990)
#' is used to calculate the expected number of events under the alternative
#' (see also Lakatos & Lan, Statistics in Medicine, 1992). These formulas are generalized
#' to piecewise survival times and non-constant piecewise accrual over time.\cr
#'
#' Optional argument \code{accountForObservationTimes}: if \code{accountForObservationTimes = FALSE},
#' only the event rates are used for the calculation of the maximum number of subjects.
#'
#' @template details_piecewise_survival
#'
#' @template details_piecewise_accrual
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family sample size functions
#'
#' @template examples_get_sample_size_survival
#'
#' @export
#'
getSampleSizeSurvival <- function(design = NULL, ...,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0 = 1, # C_THETA_H0_SURVIVAL_DEFAULT
pi1 = NA_real_,
pi2 = NA_real_,
lambda1 = NA_real_,
lambda2 = NA_real_,
median1 = NA_real_,
median2 = NA_real_,
kappa = 1,
hazardRatio = NA_real_,
piecewiseSurvivalTime = NA_real_,
allocationRatioPlanned = NA_real_, # C_ALLOCATION_RATIO_DEFAULT
eventTime = 12, # C_EVENT_TIME_DEFAULT
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
followUpTime = NA_real_,
maxNumberOfSubjects = NA_real_,
dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
dropoutTime = 12 # C_DROP_OUT_TIME_DEFAULT
) {
if (is.null(design)) {
design <- .getDefaultDesign(..., type = "sampleSize", ignore = c("accountForObservationTimes"))
.warnInCaseOfUnknownArguments(
functionName = "getSampleSizeSurvival",
ignore = c(.getDesignArgumentsToIgnoreAtUnknownArgumentCheck(
design,
powerCalculationEnabled = FALSE
), "accountForObservationTimes"), ...
)
} else {
.assertIsTrialDesign(design)
.warnInCaseOfUnknownArguments(
functionName = "getSampleSizeSurvival", ...,
ignore = c("accountForObservationTimes")
)
.warnInCaseOfTwoSidedPowerArgument(...)
}
if (!is.na(maxNumberOfSubjects) && maxNumberOfSubjects == 0) {
maxNumberOfSubjects <- NA_real_
}
# identify accrual time case
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects,
showWarnings = FALSE
)
accrualSetup$.validate()
accountForObservationTimes <- .getOptionalArgument("accountForObservationTimes", ...)
if (is.null(accountForObservationTimes)) {
accountForObservationTimes <- TRUE
}
if (!accrualSetup$maxNumberOfSubjectsCanBeCalculatedDirectly &&
accrualSetup$followUpTimeMustBeUserDefined) {
if (is.na(followUpTime)) {
if (accrualSetup$piecewiseAccrualEnabled && !accrualSetup$endOfAccrualIsUserDefined) {
stop(
C_EXCEPTION_TYPE_MISSING_ARGUMENT,
"'followUpTime', 'maxNumberOfSubjects' or end of accrual must be defined"
)
}
stop(
C_EXCEPTION_TYPE_MISSING_ARGUMENT,
"'followUpTime' or 'maxNumberOfSubjects' must be defined"
)
}
if (followUpTime == Inf) {
followUpTime <- 1e12
}
if (!any(is.na(hazardRatio)) && !is.na(thetaH0)) {
.assertIsValidHazardRatio(hazardRatio, thetaH0)
}
pwst <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda1 = lambda1, lambda2 = lambda2,
pi1 = pi1, pi2 = pi2,
median1 = median1, median2 = median2,
hazardRatio = hazardRatio, eventTime = eventTime, kappa = kappa,
.silent = TRUE
)
paramName <- NULL
if (!pwst$piecewiseSurvivalEnabled) {
if (pwst$.getParameterType("pi1") == C_PARAM_USER_DEFINED ||
pwst$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE ||
pwst$.getParameterType("pi2") == C_PARAM_USER_DEFINED) {
paramName <- "pi1"
} else if (pwst$.getParameterType("lambda1") == C_PARAM_USER_DEFINED ||
pwst$.getParameterType("lambda2") == C_PARAM_USER_DEFINED) {
paramName <- "lambda1"
} else if (pwst$.getParameterType("hazardRatio") == C_PARAM_USER_DEFINED) {
paramName <- "hazardRatio"
} else if (pwst$.getParameterType("median1") == C_PARAM_USER_DEFINED ||
pwst$.getParameterType("median2") == C_PARAM_USER_DEFINED) {
paramName <- "median1"
}
} else if (pwst$.getParameterType("hazardRatio") == C_PARAM_USER_DEFINED) {
paramName <- "hazardRatio"
}
if (!is.null(paramName)) {
paramValue <- pwst[[paramName]]
if (!is.null(paramValue) && length(paramValue) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the calculation of 'maxNumberOfSubjects' for given 'followUpTime' ",
"is only available for a single '", paramName, "'; ",
paramName, " = ", .arrayToString(
paramValue,
vectorLookAndFeelEnabled = TRUE
)
)
}
}
hr <- hazardRatio
if (all(is.na(hazardRatio))) {
hr <- pwst$hazardRatio
}
if (all(is.na(hazardRatio))) {
.assertIsValidHazardRatio(hr, thetaH0)
}
maxNumberOfSubjectsTarget <- NA_real_
withCallingHandlers(
{
# search for accrual time that provides a result
at <- accrualSetup$accrualTime
additionalAccrual <- 1
searchAccrualTimeEnabled <- TRUE
maxSearchIterations <- 50
maxNumberOfSubjectsLower <- NA_real_
maxNumberOfSubjectsLowerBefore <- 0
sampleSize <- NULL
expectionMessage <- NA_character_
while (searchAccrualTimeEnabled && maxSearchIterations >= 0 &&
(is.na(maxNumberOfSubjectsLower) ||
maxNumberOfSubjectsLower < maxNumberOfSubjectsLowerBefore ||
maxNumberOfSubjectsLower < 1e8)) {
tryCatch(
{
maxNumberOfSubjectsLowerBefore <- ifelse(is.na(maxNumberOfSubjectsLower),
0, maxNumberOfSubjectsLower
)
maxNumberOfSubjectsLower <- getAccrualTime(
accrualTime = c(at, at[length(at)] + additionalAccrual),
accrualIntensity = accrualSetup$accrualIntensity,
accrualIntensityType = accrualIntensityType
)$maxNumberOfSubjects
additionalAccrual <- 2 * additionalAccrual
sampleSize <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1, median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsLower,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2, dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)
searchAccrualTimeEnabled <- FALSE
},
error = function(e) {
expectionMessage <<- e$message
}
)
maxSearchIterations <- maxSearchIterations - 1
}
if (is.null(sampleSize) || is.na(sampleSize$followUpTime)) {
if (!is.na(expectionMessage) && grepl("'allocationRatioPlanned' > 0", expectionMessage)) {
stop(expectionMessage, call. = FALSE)
}
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"'additionalAccrual' could not be found, change accrual time specification",
call. = FALSE
)
}
# define lower bound for maxNumberOfSubjects
maxNumberOfSubjectsLower <- ceiling(max(na.omit(c(
sampleSize$eventsFixed,
as.vector(sampleSize$cumulativeEventsPerStage)
))))
if (is.na(maxNumberOfSubjectsLower)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'maxNumberOfSubjectsLower' could not be found",
call. = FALSE
)
}
# check whether accrual time already fulfills requirement
# (followUpTime < given value) or need to be increased,
# then define upper bound for maxNumberOfSubjects
maxSearchIterations <- 50
maxNumberOfSubjectsUpper <- NA_real_
fut <- sampleSize$followUpTime
iterations <- 1
while (fut <= followUpTime) {
fut <- 2 * abs(fut)
iterations <- iterations + 1
if (iterations > 50) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"search algorithm failed to end",
call. = FALSE
)
}
}
while (!is.na(fut) && fut > followUpTime && maxSearchIterations >= 0) {
maxNumberOfSubjectsUpper <- getAccrualTime(
accrualTime = c(at, at[length(at)] + additionalAccrual),
accrualIntensity = accrualSetup$accrualIntensity,
accrualIntensityType = accrualIntensityType
)$maxNumberOfSubjects
additionalAccrual <- 2 * additionalAccrual
fut <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1, median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsUpper,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2, dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)$followUpTime
maxSearchIterations <- maxSearchIterations - 1
}
if (is.na(maxNumberOfSubjectsUpper)) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"'maxNumberOfSubjectsUpper' could not be found ",
"(fut = ", fut, ", followUpTime = ", followUpTime, ")",
call. = FALSE
)
}
# use maxNumberOfSubjectsLower and maxNumberOfSubjectsUpper to find end of accrual
if (dropoutRate1 != 0 || dropoutRate2 != 0) {
# Adjust lower bound for given dropouts assuming exponential distribution
if (is.na(allocationRatioPlanned)) {
allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
}
maxNumberOfSubjectsLower <- maxNumberOfSubjectsLower /
((allocationRatioPlanned * (1 - dropoutRate1)^(
accrualSetup$accrualTime[length(accrualSetup$accrualTime)] / dropoutTime) +
(1 - dropoutRate2)^(accrualSetup$accrualTime[length(accrualSetup$accrualTime)] / dropoutTime)) /
(allocationRatioPlanned + 1))
prec <- 1
maxSearchIterations <- 50
while (prec > 1e-04 && maxSearchIterations >= 0) {
maxNumberOfSubjectsTarget <- (maxNumberOfSubjectsLower + maxNumberOfSubjectsUpper) / 2
fut <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsTarget,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = FALSE
)$followUpTime
ifelse(fut <= followUpTime,
maxNumberOfSubjectsUpper <- maxNumberOfSubjectsTarget,
maxNumberOfSubjectsLower <- maxNumberOfSubjectsTarget
)
prec <- maxNumberOfSubjectsUpper - maxNumberOfSubjectsLower
maxSearchIterations <- maxSearchIterations - 1
}
} else {
maxNumberOfSubjectsTarget <- .getOneDimensionalRootBisectionMethod(
fun = function(x) {
fut <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = x,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = FALSE
)$followUpTime
return(followUpTime - fut)
},
lower = maxNumberOfSubjectsLower,
upper = maxNumberOfSubjectsUpper,
tolerance = 1e-04,
acceptResultsOutOfTolerance = TRUE,
maxSearchIterations = 50,
direction = 0,
suppressWarnings = FALSE,
callingFunctionInformation = "getSampleSizeSurvival"
)
}
},
warning = function(w) {
invokeRestart("muffleWarning")
}
)
if (is.na(maxNumberOfSubjectsTarget)) {
stop(
C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"failed to calculate 'maxNumberOfSubjects' by given 'followUpTime' ",
"(lower = ", maxNumberOfSubjectsLower, ", upper = ", maxNumberOfSubjectsUpper, ")"
)
}
sampleSizeSurvival <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime, lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsTarget,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = FALSE
)
sampleSizeSurvival$.setParameterType("followUpTime", C_PARAM_USER_DEFINED)
sampleSizeSurvival$.accrualTime <- accrualSetup
if (!is.na(sampleSizeSurvival$followUpTime)) {
if (followUpTime == 1e12) {
followUpTime <- Inf
}
if (sampleSizeSurvival$followUpTime >= -1e-02 && sampleSizeSurvival$followUpTime <= 1e-02) {
sampleSizeSurvival$followUpTime <- 0
}
if (sampleSizeSurvival$followUpTime < followUpTime - 1e-02 ||
sampleSizeSurvival$followUpTime > followUpTime + 1e-02) {
sampleSizeSurvival$.setParameterType("followUpTime", C_PARAM_GENERATED)
warning("User defined 'followUpTime' (", followUpTime, ") ignored because ",
"follow-up time is ", round(sampleSizeSurvival$followUpTime, 4),
call. = FALSE
)
}
}
sampleSizeSurvival$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
sampleSizeSurvival$.setParameterType("accrualTime", C_PARAM_GENERATED)
return(sampleSizeSurvival)
}
return(.getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualTime,
accrualIntensity = accrualIntensity, accrualIntensityType = accrualIntensityType,
kappa = kappa, piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = followUpTime, maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio
))
}
#' @title
#' Get Power Survival
#'
#' @description
#' Returns the power, stopping probabilities, and expected sample size for testing
#' the hazard ratio in a two treatment groups survival design.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_typeOfComputation
#' @inheritParams param_allocationRatioPlanned
#' @inheritParams param_thetaH0
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_pi1_survival
#' @inheritParams param_pi2_survival
#' @inheritParams param_median1
#' @inheritParams param_median2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_directionUpper
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_eventTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param maxNumberOfEvents \code{maxNumberOfEvents > 0} is the maximum number of events, it determines
#' the power of the test and needs to be specified.
#' @inheritParams param_maxNumberOfSubjects_survival
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the power, stopping probabilities, and expected
#' sample size at given number of events and number of subjects.
#' It also calculates the time when the required events are expected under the given
#' assumptions (exponentially, piecewise exponentially, or Weibull distributed survival times
#' and constant or non-constant piecewise accrual).
#' Additionally, an allocation ratio = \code{n1 / n2} can be specified where \code{n1} and \code{n2} are the number
#' of subjects in the two treatment groups.
#'
#' The formula of Kim & Tsiatis (Biometrics, 1990)
#' is used to calculate the expected number of events under the alternative
#' (see also Lakatos & Lan, Statistics in Medicine, 1992). These formulas are generalized to piecewise survival times and
#' non-constant piecewise accrual over time.\cr
#'
#' @template details_piecewise_survival
#'
#' @template details_piecewise_accrual
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family power functions
#'
#' @template examples_get_power_survival
#'
#' @export
#'
getPowerSurvival <- function(design = NULL, ...,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0 = 1, # C_THETA_H0_SURVIVAL_DEFAULT
directionUpper = NA,
pi1 = NA_real_,
pi2 = NA_real_,
lambda1 = NA_real_,
lambda2 = NA_real_,
median1 = NA_real_,
median2 = NA_real_,
kappa = 1,
hazardRatio = NA_real_,
piecewiseSurvivalTime = NA_real_,
allocationRatioPlanned = 1, # C_ALLOCATION_RATIO_DEFAULT
eventTime = 12, # C_EVENT_TIME_DEFAULT
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
maxNumberOfSubjects = NA_real_,
maxNumberOfEvents = NA_real_,
dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
dropoutTime = 12 # C_DROP_OUT_TIME_DEFAULT
) {
if (is.null(design)) {
design <- .getDefaultDesign(..., type = "power")
.warnInCaseOfUnknownArguments(
functionName = "getPowerSurvival",
ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = TRUE), ...
)
} else {
.assertIsTrialDesign(design)
.warnInCaseOfUnknownArguments(functionName = "getPowerSurvival", ...)
.warnInCaseOfTwoSidedPowerArgument(...)
.warnInCaseOfTwoSidedPowerIsDisabled(design)
}
designPlan <- .createDesignPlanSurvival(
objectType = "power",
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0,
pi2 = pi2,
pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = TRUE,
eventTime = eventTime,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2,
lambda1 = lambda1,
median1 = median1,
median2 = median2,
directionUpper = directionUpper,
maxNumberOfEvents = maxNumberOfEvents,
maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)
.assertIsSingleNumber(allocationRatioPlanned, "allocationRatioPlanned")
.assertIsInOpenInterval(allocationRatioPlanned, "allocationRatioPlanned", 0, C_ALLOCATION_RATIO_MAXIMUM)
if (designPlan$typeOfComputation == "Schoenfeld") {
theta <- sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(log(designPlan$hazardRatio / thetaH0))
} else if (designPlan$typeOfComputation == "Freedman") {
theta <- sqrt(allocationRatioPlanned) * (designPlan$hazardRatio - 1) /
(allocationRatioPlanned * designPlan$hazardRatio + 1)
} else if (designPlan$typeOfComputation == "HsiehFreedman") {
theta <- sqrt(4 * allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(designPlan$hazardRatio - 1) / (designPlan$hazardRatio + 1)
}
if (!is.na(designPlan$directionUpper) && !designPlan$directionUpper) {
theta <- -theta
}
powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(
design = design, theta = theta, nMax = maxNumberOfEvents
)
kMax <- design$kMax
sided <- design$sided
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
numberOfResults <- length(designPlan$hazardRatio)
} else {
numberOfResults <- length(designPlan$pi1)
}
stoppingProbs <- matrix(NA_real_, kMax, numberOfResults)
designPlan$analysisTime <- matrix(NA_real_, kMax, numberOfResults)
designPlan$numberOfSubjects <- matrix(NA_real_, kMax, numberOfResults)
designPlan$studyDuration <- rep(NA_real_, numberOfResults)
designPlan$expectedNumberOfSubjects <- rep(NA_real_, numberOfResults)
cumulativeEventsPerStage <- maxNumberOfEvents * design$informationRates
parameterNames <- c(
"analysisTime",
"numberOfSubjects",
"studyDuration",
"expectedNumberOfSubjects",
"cumulativeEventsPerStage"
)
phi <- -c(log(1 - dropoutRate1), log(1 - dropoutRate2)) / dropoutTime
lambda1 <- designPlan$lambda1
if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
lambda1 <- rep(NA_real_, numberOfResults)
}
for (i in 1:numberOfResults) {
# Analysis times
up <- 2
iterate <- 1
while (cumulativeEventsPerStage[kMax] / designPlan$maxNumberOfSubjects > .getEventProbabilities(
time = up, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i], phi = phi,
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)) {
up <- 2 * up
iterate <- iterate + 1
if (iterate > 50) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ") ",
"is too small to reach maximum number of events ",
"(presumably due to drop-out rates)"
)
}
}
for (j in 1:kMax) {
designPlan$analysisTime[j, i] <- .getOneDimensionalRoot(
function(x) {
cumulativeEventsPerStage[j] / designPlan$maxNumberOfSubjects -
.getEventProbabilities(
time = x,
accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i], phi = phi,
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = up, tolerance = 1e-06,
callingFunctionInformation = "getPowerSurvival"
)
if (is.na(designPlan$analysisTime[j, i])) {
warning("Cannot calculate analysis time at stage ", j, ": ",
"'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ") is too ",
"small to reach maximum number of events",
call. = FALSE
)
}
}
if (kMax > 1) {
designPlan$numberOfSubjects[, i] <- .getNumberOfSubjects(
time = designPlan$analysisTime[, i],
accrualTime = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = designPlan$maxNumberOfSubjects
)
powerAndAverageSampleNumber$futilityPerStage[is.na(
powerAndAverageSampleNumber$futilityPerStage[, i]
), i] <- 0
stoppingProbs[, i] <- powerAndAverageSampleNumber$rejectPerStage[, i] +
c(powerAndAverageSampleNumber$futilityPerStage[, i], 0)
stoppingProbs[kMax, i] <- 1 - sum(stoppingProbs[1:(kMax - 1), i])
designPlan$studyDuration[i] <- designPlan$analysisTime[, i] %*% stoppingProbs[, i]
designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
designPlan$expectedNumberOfSubjects[i] <-
designPlan$numberOfSubjects[, i] %*% stoppingProbs[, i]
}
}
if (kMax == 1) {
designPlan$expectedNumberOfSubjects <- .getNumberOfSubjects(
time = designPlan$analysisTime[1, ],
accrualTime = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = designPlan$maxNumberOfSubjects
)
designPlan$numberOfSubjects <- matrix(designPlan$expectedNumberOfSubjects, nrow = 1)
}
designPlan$cumulativeEventsPerStage <- matrix(cumulativeEventsPerStage, ncol = 1)
designPlan$.setParameterType("cumulativeEventsPerStage", C_PARAM_GENERATED)
.addDeprecatedFieldValues(designPlan, "eventsPerStage", designPlan$cumulativeEventsPerStage)
designPlan$expectedNumberOfEvents <- powerAndAverageSampleNumber$averageSampleNumber
designPlan$overallReject <- powerAndAverageSampleNumber$overallReject
designPlan$rejectPerStage <- powerAndAverageSampleNumber$rejectPerStage
designPlan$futilityStop <- powerAndAverageSampleNumber$overallFutility
designPlan$futilityPerStage <- powerAndAverageSampleNumber$futilityPerStage
designPlan$earlyStop <- powerAndAverageSampleNumber$overallEarlyStop
parameterNames <- c(
parameterNames,
"expectedNumberOfEvents",
"overallReject",
"rejectPerStage",
"futilityStop",
"futilityPerStage",
"earlyStop"
)
for (parameterName in parameterNames) {
designPlan$.setParameterType(parameterName, C_PARAM_GENERATED)
}
if (kMax == 1L) {
designPlan$.setParameterType("numberOfSubjects", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("cumulativeEventsPerStage", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("futilityStop", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("earlyStop", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("rejectPerStage", C_PARAM_NOT_APPLICABLE)
}
if (!any(is.na(designPlan$analysisTime)) && !any(is.na(designPlan$accrualTime))) {
designPlan$followUpTime <- designPlan$analysisTime[kMax, ] -
designPlan$accrualTime[length(designPlan$accrualTime)]
designPlan$.setParameterType("followUpTime", C_PARAM_GENERATED)
}
designPlan$singleEventsPerStage <- .getSingleEventsPerStage(designPlan$cumulativeEventsPerStage)
.addEffectScaleBoundaryDataToDesignPlan(designPlan)
.addStudyDurationToDesignPlan(designPlan)
.hideFutilityStopsIfNotApplicable(designPlan)
return(designPlan)
}
Try the rpact package in your browser
Any scripts or data that you put into this service are public.
rpact documentation built on May 29, 2024, 11:20 a.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 getPowerSurvival getSampleSizeSurvival getNumberOfSubjects getEventProbabilities .getNumberOfSubjects .getSampleSizeSequentialSurvival .getSampleSizeFixedSurvival .initDesignPlanSurvival .initDesignPlanSurvivalByPiecewiseSurvivalTimeObject .createDesignPlanSurvival .getSampleSizeSurvival .calculateSampleSizeSurvival .getEffectScaleBoundaryDataSurvival .getEventProbabilities .getEventsFixed .isUserDefinedMaxNumberOfSubjects .getEventProbabilitiesOverall .getEventProbabilitiesGroupwise .getEventProbabilityFunctionVec .getEventProbabilityFunction .getLambda .addStudyDurationToDesignPlan
Documented in getEventProbabilities getNumberOfSubjects getPowerSurvival getSampleSizeSurvival
## |
## | *Sample size and power survival*
## |
## | 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: 7662 $
## | Last changed: $Date: 2024-02-23 12:42:26 +0100 (Fr, 23 Feb 2024) $
## | Last changed by: $Author: pahlke $
## |
#' @include f_core_utilities.R
NULL
.addStudyDurationToDesignPlan <- function(designPlan) {
if (!designPlan$accountForObservationTimes) {
return(invisible())
}
kMax <- designPlan$.design$kMax
if (kMax == 1) {
designPlan$studyDuration <- designPlan$analysisTime[1, ]
designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
designPlan$maxStudyDuration <- designPlan$studyDuration
} else {
designPlan$maxStudyDuration <- designPlan$analysisTime[kMax, ]
designPlan$.setParameterType("maxStudyDuration", C_PARAM_GENERATED)
}
}
.getLambda <- function(..., groupNumber, lambda2, lambda1, hazardRatio, kappa) {
if (groupNumber == 1) {
if (!any(is.na(lambda1))) {
return(lambda1)
}
lambda2 <- lambda2 * hazardRatio^(1 / kappa)
}
return(lambda2)
}
.getEventProbabilityFunction <- function(..., time, piecewiseLambda, piecewiseSurvivalTime, phi, kappa) {
if (length(piecewiseLambda) == 1) {
if (kappa != 1 && phi > 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "Weibull distribution cannot ",
"be used together with specified dropout rate (use simulation instead)",
call. = FALSE
)
}
return(piecewiseLambda / (piecewiseLambda + phi) *
pweibull(time, shape = kappa, scale = 1 / (piecewiseLambda + phi), lower.tail = TRUE, log.p = FALSE))
}
if (length(piecewiseSurvivalTime) != length(piecewiseLambda)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of 'piecewiseSurvivalTime' (", .arrayToString(piecewiseSurvivalTime),
") must be equal to length of 'piecewiseLambda' (", .arrayToString(piecewiseLambda), ")"
)
}
piecewiseSurvivalTime <- .getPiecewiseExpStartTimesWithoutLeadingZero(piecewiseSurvivalTime)
if (kappa != 1) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"Weibull distribution cannot be used for piecewise survival definition",
call. = FALSE
)
}
len <- length(piecewiseSurvivalTime)
for (i in 1:len) {
if (i == 1) {
if (time <= piecewiseSurvivalTime[1]) {
return(piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
(1 - exp(-((piecewiseLambda[1] + phi) * time))))
}
} else if (i == 2) {
cdfPart <- piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
(1 - exp(-((piecewiseLambda[1] + phi) * piecewiseSurvivalTime[1])))
if (time <= piecewiseSurvivalTime[2]) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
cdf <- cdfPart + piecewiseLambda[2] / (piecewiseLambda[2] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[1]) - exp(-piecewiseLambda[2] *
(time - piecewiseSurvivalTime[1]) - phi * time))
return(cdf)
}
} else if (i == 3) {
cdfPart <- cdfPart + piecewiseLambda[2] / (piecewiseLambda[2] + phi) *
exp(-piecewiseLambda[1] * piecewiseSurvivalTime[1]) * (
exp(-phi * piecewiseSurvivalTime[1]) - exp(-piecewiseLambda[2] *
(piecewiseSurvivalTime[2] - piecewiseSurvivalTime[1]) - phi * piecewiseSurvivalTime[2]))
if (time <= piecewiseSurvivalTime[3]) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
piecewiseLambda[2] * (piecewiseSurvivalTime[2] - piecewiseSurvivalTime[1])
cdf <- cdfPart + piecewiseLambda[3] / (piecewiseLambda[3] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[2]) - exp(-piecewiseLambda[3] *
(time - piecewiseSurvivalTime[2]) - phi * time))
return(cdf)
}
} else if (i > 3) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
sum(piecewiseLambda[2:(i - 2)] * (piecewiseSurvivalTime[2:(i - 2)] -
piecewiseSurvivalTime[1:(i - 3)]))
cdfPart <- cdfPart + piecewiseLambda[i - 1] / (piecewiseLambda[i - 1] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[i - 2]) - exp(-piecewiseLambda[i - 1] *
(piecewiseSurvivalTime[i - 1] - piecewiseSurvivalTime[i - 2]) - phi * piecewiseSurvivalTime[i - 1]))
if (time <= piecewiseSurvivalTime[i]) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
sum(piecewiseLambda[2:(i - 1)] * (piecewiseSurvivalTime[2:(i - 1)] - piecewiseSurvivalTime[1:(i - 2)]))
cdf <- cdfPart + piecewiseLambda[i] / (piecewiseLambda[i] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[i - 1]) - exp(-piecewiseLambda[i] *
(time - piecewiseSurvivalTime[i - 1]) - phi * time))
return(cdf)
}
}
}
if (len == 1) {
cdfPart <- piecewiseLambda[1] / (piecewiseLambda[1] + phi) *
(1 - exp(-((piecewiseLambda[1] + phi) * piecewiseSurvivalTime[1])))
} else if (len == 2) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
cdfPart <- cdfPart + piecewiseLambda[len] / (piecewiseLambda[len] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[len - 1]) - exp(-piecewiseLambda[len] *
(piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1]) - phi * piecewiseSurvivalTime[len]))
} else {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1] +
sum(piecewiseLambda[2:(len - 1)] * (piecewiseSurvivalTime[2:(len - 1)] - piecewiseSurvivalTime[1:(len - 2)]))
cdfPart <- cdfPart + piecewiseLambda[len] / (piecewiseLambda[len] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[len - 1]) - exp(-piecewiseLambda[len] *
(piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1]) - phi * piecewiseSurvivalTime[len]))
}
if (len == 1) {
cdfFactor <- piecewiseLambda[1] * piecewiseSurvivalTime[1]
} else {
cdfFactor <- cdfFactor + piecewiseLambda[len] * (piecewiseSurvivalTime[len] - piecewiseSurvivalTime[len - 1])
}
cdf <- cdfPart + piecewiseLambda[len + 1] / (piecewiseLambda[len + 1] + phi) * exp(-cdfFactor) * (
exp(-phi * piecewiseSurvivalTime[len]) - exp(-piecewiseLambda[len + 1] *
(time - piecewiseSurvivalTime[len]) - phi * time))
return(cdf)
}
.getEventProbabilityFunctionVec <- function(..., timeVector, piecewiseLambda, piecewiseSurvivalTime, phi, kappa) {
result <- c()
for (time in timeVector) {
result <- c(result, .getEventProbabilityFunction(
time = time, piecewiseLambda = piecewiseLambda,
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa
))
}
return(result)
}
.getEventProbabilitiesGroupwise <- function(
..., time, accrualTimeVector, accrualIntensity, lambda2,
lambda1, piecewiseSurvivalTime, phi, kappa, allocationRatioPlanned, hazardRatio) {
.assertIsSingleNumber(time, "time")
if (length(accrualTimeVector) > 1 && accrualTimeVector[1] == 0) {
accrualTimeVector <- accrualTimeVector[2:length(accrualTimeVector)]
}
accrualTimeVectorLength <- length(accrualTimeVector)
densityIntervals <- accrualTimeVector
if (accrualTimeVectorLength > 1) {
densityIntervals[2:accrualTimeVectorLength] <-
accrualTimeVector[2:accrualTimeVectorLength] -
accrualTimeVector[1:(accrualTimeVectorLength - 1)]
}
if (length(densityIntervals) > 1 && length(accrualIntensity) > 1 &&
length(densityIntervals) != length(accrualIntensity)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'densityIntervals' (", .arrayToString(densityIntervals),
") and 'accrualIntensity' (", .arrayToString(accrualIntensity), ") must have same length"
)
}
densityVector <- accrualIntensity / sum(densityIntervals * accrualIntensity)
eventProbs <- rep(NA_real_, 2)
for (k in 1:accrualTimeVectorLength) {
if (time <= accrualTimeVector[k]) {
for (groupNumber in c(1, 2)) { # two groups: 1 = treatment, 2 = control
lambdaTemp <- .getLambda(
groupNumber = groupNumber,
lambda2 = lambda2,
lambda1 = lambda1,
hazardRatio = hazardRatio,
kappa = kappa
)
inner <- function(x) {
.getEventProbabilityFunctionVec(
timeVector = x, piecewiseLambda = lambdaTemp,
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi[groupNumber], kappa = kappa
)
}
timeValue1 <- 0
if (k > 1) {
timeValue1 <- time - accrualTimeVector[1]
}
eventProbs[groupNumber] <- densityVector[1] * stats::integrate(inner, timeValue1, time)$value
if (k > 2) {
for (j in 2:(k - 1)) {
eventProbs[groupNumber] <- eventProbs[groupNumber] +
densityVector[j] * stats::integrate(
inner, time - accrualTimeVector[j],
time - accrualTimeVector[j - 1]
)$value
}
}
if (k > 1) {
eventProbs[groupNumber] <- eventProbs[groupNumber] +
densityVector[k] * stats::integrate(inner, 0, time - accrualTimeVector[k - 1])$value
}
}
return(eventProbs)
}
}
for (groupNumber in c(1, 2)) {
lambdaTemp <- .getLambda(
groupNumber = groupNumber,
lambda2 = lambda2,
lambda1 = lambda1,
hazardRatio = hazardRatio,
kappa = kappa
)
inner <- function(x) {
.getEventProbabilityFunctionVec(
timeVector = x, piecewiseLambda = lambdaTemp,
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi[groupNumber], kappa = kappa
)
}
eventProbs[groupNumber] <- densityVector[1] *
stats::integrate(inner, time - accrualTimeVector[1], time)$value
if (accrualTimeVectorLength > 1) {
for (j in (2:accrualTimeVectorLength)) {
eventProbs[groupNumber] <- eventProbs[groupNumber] +
densityVector[j] * stats::integrate(
inner, time - accrualTimeVector[j],
time - accrualTimeVector[j - 1]
)$value
}
}
}
return(eventProbs)
}
.getEventProbabilitiesOverall <- function(eventProbs, allocationRatioPlanned) {
return((allocationRatioPlanned * eventProbs[1] + eventProbs[2]) / (1 + allocationRatioPlanned))
}
.isUserDefinedMaxNumberOfSubjects <- function(designPlan) {
if (!is.null(designPlan) && length(designPlan$.getParameterType("maxNumberOfSubjects")) > 0) {
if (designPlan$.getParameterType("maxNumberOfSubjects") == C_PARAM_USER_DEFINED) {
return(TRUE)
}
}
return(!is.null(designPlan$maxNumberOfSubjects) && length(designPlan$maxNumberOfSubjects) == 1 &&
!is.na(designPlan$maxNumberOfSubjects) && designPlan$maxNumberOfSubjects > 0)
}
.getEventsFixed <- function(..., typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
twoSidedPower, alpha, beta, sided, hazardRatio, thetaH0, allocationRatioPlanned) {
typeOfComputation <- match.arg(typeOfComputation)
if (typeOfComputation == "Schoenfeld") {
eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 /
(log(hazardRatio) - log(thetaH0))^2 *
(1 + allocationRatioPlanned)^2 / allocationRatioPlanned
if (twoSidedPower && (sided == 2)) {
up <- 2 * eventsFixed
eventsFixed <- .getOneDimensionalRoot(
function(n) {
return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
(log(hazardRatio) - log(thetaH0)) * sqrt(allocationRatioPlanned) /
(1 + allocationRatioPlanned)) -
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
(log(hazardRatio) - log(thetaH0)) * sqrt(allocationRatioPlanned) /
(1 + allocationRatioPlanned)) - beta)
},
lower = 0.001, upper = up, tolerance = 1e-04,
callingFunctionInformation = ".getEventsFixed"
)
}
return(eventsFixed)
}
if (typeOfComputation == "Freedman") {
eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 *
(1 + hazardRatio * allocationRatioPlanned)^2 / (1 - hazardRatio)^2 /
allocationRatioPlanned
if (twoSidedPower && (sided == 2)) {
up <- 2 * eventsFixed
eventsFixed <- .getOneDimensionalRoot(
function(n) {
return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
sqrt(allocationRatioPlanned) * (1 - hazardRatio) /
(1 + allocationRatioPlanned * hazardRatio)) -
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
sqrt(allocationRatioPlanned) * (1 - hazardRatio) /
(1 + allocationRatioPlanned * hazardRatio)) - beta)
},
lower = 0.001, upper = up, tolerance = 1e-04,
callingFunctionInformation = ".getEventsFixed"
)
}
return(eventsFixed)
}
if (typeOfComputation == "HsiehFreedman") {
eventsFixed <- (.getOneMinusQNorm(alpha / sided) + .getOneMinusQNorm(beta))^2 *
(1 + hazardRatio)^2 / (1 - hazardRatio)^2 *
(1 + allocationRatioPlanned)^2 / (4 * allocationRatioPlanned)
if (twoSidedPower && sided == 2) {
up <- 2 * eventsFixed
eventsFixed <- .getOneDimensionalRoot(
function(n) {
return(stats::pnorm(.getOneMinusQNorm(alpha / 2) - sqrt(n) *
2 * sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(1 - hazardRatio) / (1 + hazardRatio)) -
stats::pnorm(-.getOneMinusQNorm(alpha / 2) - sqrt(n) *
2 * sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(1 - hazardRatio) / (1 + hazardRatio)) - beta)
},
lower = 0.001, upper = up, tolerance = 1e-04,
callingFunctionInformation = ".getEventsFixed"
)
}
return(eventsFixed)
}
}
.getEventProbabilities <- function(..., time, accrualTimeVector, accrualIntensity, lambda2,
lambda1, piecewiseSurvivalTime, phi, kappa, allocationRatioPlanned, hazardRatio) {
eventProbs <- .getEventProbabilitiesGroupwise(
time = time, accrualTimeVector = accrualTimeVector,
accrualIntensity = accrualIntensity, lambda2 = lambda2,
lambda1 = lambda1, piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio
)
return(.getEventProbabilitiesOverall(eventProbs, allocationRatioPlanned))
}
.getEffectScaleBoundaryDataSurvival <- function(designPlan) {
design <- designPlan$.design
thetaH0 <- designPlan$thetaH0
cumulativeEventsPerStage <- designPlan$cumulativeEventsPerStage
allocationRatioPlanned <- designPlan$allocationRatioPlanned
directionUpper <- designPlan$directionUpper
if (design$kMax == 1) {
nParameters <- length(cumulativeEventsPerStage)
} else {
nParameters <- ncol(cumulativeEventsPerStage)
}
directionUpper[is.na(directionUpper)] <- TRUE
if (length(allocationRatioPlanned) == 1) {
allocationRatioPlanned <- rep(allocationRatioPlanned, nParameters)
}
futilityBounds <- design$futilityBounds
futilityBounds[!is.na(futilityBounds) & futilityBounds <= C_FUTILITY_BOUNDS_DEFAULT] <- NA_real_
criticalValues <- design$criticalValues
criticalValuesEffectScaleUpper <- matrix(, nrow = design$kMax, ncol = nParameters)
criticalValuesEffectScaleLower <- matrix(, nrow = design$kMax, ncol = nParameters)
futilityBoundsEffectScaleUpper <- matrix(, nrow = design$kMax - 1, ncol = nParameters)
futilityBoundsEffectScaleLower <- matrix(, nrow = design$kMax - 1, ncol = nParameters)
for (j in (1:nParameters)) {
if (design$sided == 1) {
criticalValuesEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * criticalValues *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[, j])))
} else {
criticalValuesEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * criticalValues *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[, j])))
criticalValuesEffectScaleLower[, j] <- thetaH0 * (exp(-(2 * directionUpper[j] - 1) * criticalValues *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[, j])))
}
if (!.isTrialDesignFisher(design) && !all(is.na(futilityBounds))) {
futilityBoundsEffectScaleUpper[, j] <- thetaH0 * (exp((2 * directionUpper[j] - 1) * futilityBounds *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[1:(design$kMax - 1), j])))
}
if (!.isTrialDesignFisher(design) && design$sided == 2 && design$kMax > 1 &&
(design$typeOfDesign == C_TYPE_OF_DESIGN_PT || !is.null(design$typeBetaSpending) && design$typeBetaSpending != "none")) {
futilityBoundsEffectScaleLower[, j] <- thetaH0 * (exp(-(2 * directionUpper[j] - 1) * futilityBounds *
(1 + allocationRatioPlanned[j]) / sqrt(allocationRatioPlanned[j] *
cumulativeEventsPerStage[1:(design$kMax - 1), j])))
}
}
return(list(
criticalValuesEffectScaleUpper = matrix(criticalValuesEffectScaleUpper, nrow = design$kMax),
criticalValuesEffectScaleLower = matrix(criticalValuesEffectScaleLower, nrow = design$kMax),
futilityBoundsEffectScaleUpper = matrix(futilityBoundsEffectScaleUpper, nrow = design$kMax - 1),
futilityBoundsEffectScaleLower = matrix(futilityBoundsEffectScaleLower, nrow = design$kMax - 1)
))
}
.calculateSampleSizeSurvival <- function(designPlan) {
# Fixed
designPlan <- .getSampleSizeFixedSurvival(designPlan)
# Sequential
if (designPlan$.design$kMax > 1) {
designCharacteristics <- getDesignCharacteristics(designPlan$.design)
designPlan <- .getSampleSizeSequentialSurvival(designPlan, designCharacteristics)
}
if (designPlan$accountForObservationTimes && !any(is.na(designPlan$followUpTime)) &&
all(designPlan$followUpTime == C_FOLLOW_UP_TIME_DEFAULT)) {
designPlan$.setParameterType("followUpTime", C_PARAM_DEFAULT_VALUE)
}
.addEffectScaleBoundaryDataToDesignPlan(designPlan)
designPlan$singleEventsPerStage <- .getSingleEventsPerStage(designPlan$cumulativeEventsPerStage)
.addDeprecatedFieldValues(designPlan, "eventsPerStage", designPlan$cumulativeEventsPerStage)
paramName <- ifelse(designPlan$.design$kMax == 1, "nFixed", "maxNumberOfSubjects")
if (designPlan$.getParameterType(paramName) == C_PARAM_GENERATED &&
designPlan$.accrualTime$.getParameterType("maxNumberOfSubjects") != C_PARAM_GENERATED &&
all(designPlan$accrualIntensity < 1)) {
numberOfDefinedAccrualIntensities <- length(designPlan$accrualIntensity)
accrualTime <- designPlan$accrualTime
if (length(accrualTime) > 0 && accrualTime[1] != 0) {
accrualTime <- c(0, accrualTime)
}
if (any(designPlan$accrualIntensity < 0)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'accrualIntensityRelative' (",
.arrayToString(designPlan$accrualIntensity), ") must be >= 0"
)
}
designPlan$accrualIntensityRelative <- designPlan$accrualIntensity
if (identical(designPlan$accrualIntensityRelative, C_ACCRUAL_INTENSITY_DEFAULT)) {
designPlan$.setParameterType("accrualIntensityRelative", C_PARAM_NOT_APPLICABLE)
} else {
designPlan$.setParameterType(
"accrualIntensityRelative",
designPlan$.getParameterType("accrualIntensity")
)
}
accrualIntensityAbsolute <- numeric()
for (maxNumberOfSubjects in designPlan[[paramName]]) {
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = designPlan$accrualIntensityRelative,
accrualIntensityType = "relative",
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualIntensityAbsolute <- c(accrualIntensityAbsolute, accrualSetup$accrualIntensity)
}
designPlan$accrualIntensity <- accrualIntensityAbsolute
designPlan$.setParameterType("accrualIntensity", C_PARAM_GENERATED)
if (numberOfDefinedAccrualIntensities > 1) {
paramName <- NULL
if (designPlan$.getParameterType("pi1") == C_PARAM_USER_DEFINED ||
designPlan$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE ||
designPlan$.getParameterType("pi2") == C_PARAM_USER_DEFINED) {
paramName <- "pi1"
} else if (designPlan$.getParameterType("median1") == C_PARAM_USER_DEFINED ||
designPlan$.getParameterType("median2") == C_PARAM_USER_DEFINED) {
paramName <- "median1"
}
if (!is.null(paramName)) {
paramValue <- designPlan[[paramName]]
if (!is.null(paramValue) && length(paramValue) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the definition of relative accrual intensities ",
"(all 'accrualIntensity' values < 1) ",
"is only available for a single value ",
"(", paramName, " = ", .arrayToString(
paramValue,
vectorLookAndFeelEnabled = TRUE
), ")"
)
}
}
}
}
designPlan$maxNumberOfEvents <- designPlan$cumulativeEventsPerStage[designPlan$.design$kMax, ]
designPlan$.setParameterType("maxNumberOfEvents", C_PARAM_GENERATED)
if (!any(is.na(designPlan$followUpTime))) {
if (any(designPlan$followUpTime < -1e-02)) {
warning("Accrual duration longer than maximal study ",
"duration (time to maximal number of events); followUpTime = ",
.arrayToString(designPlan$followUpTime),
call. = FALSE
)
}
} else {
warning("Follow-up time could not be calculated for hazardRatio = ",
.arrayToString(designPlan$hazardRatio[indices]),
call. = FALSE
)
}
if (designPlan$.getParameterType("accountForObservationTimes") != C_PARAM_USER_DEFINED) {
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_NOT_APPLICABLE)
}
designPlan$.setParameterType("chi", C_PARAM_NOT_APPLICABLE)
.addStudyDurationToDesignPlan(designPlan)
return(designPlan)
}
# Hidden parameter:
# @param accountForObservationTimes If \code{accountForObservationTimes = TRUE}, the number of
# subjects is calculated assuming specific accrual and follow-up time, default is \code{TRUE}
# (see details).
# If \code{accountForObservationTimes = FALSE}, only the event rates are used for the calculation
# of the maximum number of subjects.
# \code{accountForObservationTimes} can be selected as \code{FALSE}. In this case,
# the number of subjects is calculated from the event probabilities only.
# This kind of computation does not account for the specific accrual pattern and survival distribution.
.getSampleSizeSurvival <- function(...,
design = NULL,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0 = 1,
pi2 = NA_real_,
pi1 = NA_real_,
allocationRatioPlanned = NA_real_,
accountForObservationTimes = TRUE,
eventTime = C_EVENT_TIME_DEFAULT,
accrualTime = C_ACCRUAL_TIME_DEFAULT,
accrualIntensity = C_ACCRUAL_INTENSITY_DEFAULT,
accrualIntensityType = c("auto", "absolute", "relative"),
kappa = 1,
piecewiseSurvivalTime = NA_real_,
lambda2 = NA_real_,
lambda1 = NA_real_,
median1 = NA_real_,
median2 = NA_real_,
followUpTime = NA_real_,
maxNumberOfSubjects = NA_real_,
dropoutRate1 = 0,
dropoutRate2 = dropoutRate1,
dropoutTime = NA_real_,
hazardRatio = NA_real_,
userFunctionCallEnabled = TRUE) {
designPlan <- .createDesignPlanSurvival(
objectType = "sampleSize",
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0,
pi2 = pi2,
pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2,
lambda1 = lambda1,
median1 = median1,
median2 = median2,
followUpTime = followUpTime,
maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = userFunctionCallEnabled
)
return(.calculateSampleSizeSurvival(designPlan))
}
.createDesignPlanSurvival <- function(..., objectType = c("sampleSize", "power"),
design,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0, pi2, pi1,
allocationRatioPlanned,
accountForObservationTimes,
eventTime,
accrualTime,
accrualIntensity,
accrualIntensityType,
kappa,
piecewiseSurvivalTime,
lambda2,
lambda1,
median1,
median2,
followUpTime = NA_real_,
directionUpper = NA,
maxNumberOfEvents = NA_real_,
maxNumberOfSubjects,
dropoutRate1,
dropoutRate2,
dropoutTime,
hazardRatio,
userFunctionCallEnabled = TRUE) {
objectType <- match.arg(objectType)
typeOfComputation <- .matchArgument(typeOfComputation, "Schoenfeld")
if (userFunctionCallEnabled) {
.assertIsTrialDesignInverseNormalOrGroupSequential(design)
.assertIsValidAlphaAndBeta(design$alpha, design$beta)
.assertIsValidSidedParameter(design$sided)
.assertIsSingleLogical(accountForObservationTimes, "accountForObservationTimes", naAllowed = TRUE)
.assertIsSingleNumber(thetaH0, "thetaH0")
.assertIsValidThetaH0(thetaH0, endpoint = "survival", groups = 2)
.assertIsValidKappa(kappa)
directionUpper <- .assertIsValidDirectionUpper(directionUpper,
design$sided, objectType,
userFunctionCallEnabled = TRUE
)
if (objectType == "power") {
.assertIsSingleNumber(maxNumberOfEvents, "maxNumberOfEvents")
.assertIsInClosedInterval(maxNumberOfEvents, "maxNumberOfEvents",
lower = 1, upper = maxNumberOfSubjects
)
}
if (!any(is.na(pi1)) && (any(pi1 <= 0) || any(pi1 >= 1))) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"event rate 'pi1' (", .arrayToString(pi1), ") is out of bounds (0; 1)"
)
}
if (!any(is.na(pi2)) && (any(pi2 <= 0) || any(pi2 >= 1))) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"event rate 'pi2' (", .arrayToString(pi2), ") is out of bounds (0; 1)"
)
}
if (design$sided == 2 && thetaH0 != 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"two-sided case is implemented for superiority testing only (i.e., thetaH0 = 1)"
)
}
if (thetaH0 <= 0) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"null hypothesis hazard ratio is not allowed be negative or zero"
)
}
if (!(typeOfComputation %in% c("Schoenfeld", "Freedman", "HsiehFreedman"))) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"computation type ('", typeOfComputation, "') must be one of the following: ",
"'Schoenfeld', 'Freedman', or 'HsiehFreedman' "
)
}
if (typeOfComputation != "Schoenfeld" && thetaH0 != 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"Freedman test calculation is possible only for superiority testing (thetaH0 != 1)"
)
}
}
if (is.numeric(accrualTime) && all(is.na(accrualTime))) {
accrualTime <- C_ACCRUAL_TIME_DEFAULT
}
if (all(is.na(accrualIntensity))) {
accrualIntensity <- C_ACCRUAL_INTENSITY_DEFAULT
}
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualSetup$.validate()
if (is.na(allocationRatioPlanned)) {
allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
}
if (userFunctionCallEnabled) {
.assertIsValidAllocationRatioPlannedSampleSize(allocationRatioPlanned, maxNumberOfSubjects)
}
designPlan <- TrialDesignPlanSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime,
accrualTime = accrualSetup$.getAccrualTimeWithoutLeadingZero(),
accrualIntensity = accrualSetup$accrualIntensity,
kappa = kappa,
followUpTime = followUpTime,
maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)
.setValueAndParameterType(
designPlan, "allocationRatioPlanned",
allocationRatioPlanned, C_ALLOCATION_RATIO_DEFAULT
)
.setValueAndParameterType(designPlan, "dropoutRate1", dropoutRate1, C_DROP_OUT_RATE_1_DEFAULT)
.setValueAndParameterType(designPlan, "dropoutRate2", dropoutRate2, C_DROP_OUT_RATE_2_DEFAULT)
.setValueAndParameterType(designPlan, "dropoutTime", dropoutTime, C_DROP_OUT_TIME_DEFAULT)
.setValueAndParameterType(designPlan, "kappa", kappa, 1)
designPlan$.setObjectType(objectType)
designPlan$criticalValuesPValueScale <- matrix(design$stageLevels, ncol = 1)
if (design$sided == 2) {
designPlan$criticalValuesPValueScale <- designPlan$criticalValuesPValueScale * 2
designPlan$.setParameterType("criticalValuesPValueScale", C_PARAM_GENERATED)
}
if (any(design$futilityBounds > C_FUTILITY_BOUNDS_DEFAULT)) {
designPlan$futilityBoundsPValueScale <- matrix(1 - stats::pnorm(design$futilityBounds), ncol = 1)
designPlan$.setParameterType("futilityBoundsPValueScale", C_PARAM_GENERATED)
}
designPlan$.accrualTime <- accrualSetup
designPlan$totalAccrualTime <- accrualSetup$accrualTime[length(accrualSetup$accrualTime)]
if (length(accrualSetup$accrualTime) > 2) {
designPlan$.setParameterType("totalAccrualTime", C_PARAM_GENERATED)
} else {
designPlan$.setParameterType("totalAccrualTime", C_PARAM_NOT_APPLICABLE)
}
if (is.na(maxNumberOfSubjects)) {
if (!is.na(accrualSetup$maxNumberOfSubjects)) {
designPlan$maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects
designPlan$.setParameterType(
"maxNumberOfSubjects",
accrualSetup$.getParameterType("maxNumberOfSubjects")
)
}
} else if (maxNumberOfSubjects == 0) {
designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
} else {
designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_USER_DEFINED)
}
if (identical(as.integer(accrualSetup$accrualTime), C_ACCRUAL_TIME_DEFAULT) ||
identical(
as.integer(c(0L, accrualSetup$.getAccrualTimeWithoutLeadingZero())),
C_ACCRUAL_TIME_DEFAULT
)) {
designPlan$.setParameterType("accrualTime", C_PARAM_DEFAULT_VALUE)
} else {
designPlan$.setParameterType("accrualTime", accrualSetup$.getParameterType("accrualTime"))
}
if (length(designPlan$accrualIntensity) == 1 &&
designPlan$accrualIntensity == C_ACCRUAL_INTENSITY_DEFAULT) {
designPlan$.setParameterType("accrualIntensity", C_PARAM_DEFAULT_VALUE)
} else {
designPlan$.setParameterType(
"accrualIntensity",
accrualSetup$.getParameterType("accrualIntensity")
)
}
.assertIsSingleNumber(designPlan$eventTime, "eventTime")
.assertIsSingleNumber(designPlan$allocationRatioPlanned, "allocationRatioPlanned")
.assertIsSingleNumber(designPlan$kappa, "kappa")
if (objectType == "power") {
.assertIsValidMaxNumberOfSubjects(designPlan$maxNumberOfSubjects)
}
.assertIsSingleNumber(designPlan$dropoutRate1, "dropoutRate1")
.assertIsSingleNumber(designPlan$dropoutRate2, "dropoutRate2")
.assertIsSingleNumber(designPlan$dropoutTime, "dropoutTime")
if (objectType == "power") {
pi1Default <- C_PI_1_DEFAULT
} else {
pi1Default <- C_PI_1_SAMPLE_SIZE_DEFAULT
}
designPlan$.piecewiseSurvivalTime <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime, lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
hazardRatio = hazardRatio, pi1 = pi1, pi2 = pi2, eventTime = eventTime, kappa = kappa,
.pi1Default = pi1Default
)
designPlan$.setParameterType("kappa", designPlan$.piecewiseSurvivalTime$.getParameterType("kappa"))
if (designPlan$.piecewiseSurvivalTime$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE &&
length(designPlan$.piecewiseSurvivalTime$pi1) > 1 &&
length(accrualSetup$accrualIntensity) > 1 && all(accrualSetup$accrualIntensity < 1)) {
designPlan$.piecewiseSurvivalTime$pi1 <- designPlan$.piecewiseSurvivalTime$pi1[1]
warning("Only the first default 'pi1' (", designPlan$.piecewiseSurvivalTime$pi1, ") was used ",
"because the accrual intensities (", .arrayToString(accrualSetup$accrualIntensity), ") ",
"were defined relative (all accrual intensities are < 1)",
call. = FALSE
)
}
.initDesignPlanSurvival(designPlan)
designPlan$.setParameterType("followUpTime", C_PARAM_NOT_APPLICABLE)
if (designPlan$accountForObservationTimes) {
.assertIsSingleNumber(dropoutRate1, "dropoutRate1")
.assertIsSingleNumber(dropoutRate2, "dropoutRate2")
.assertIsSingleNumber(dropoutTime, "dropoutTime")
if (!is.na(dropoutTime) && dropoutTime <= 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'dropoutTime' (", dropoutTime, ") must be > 0")
}
if (dropoutRate1 < 0 || dropoutRate1 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate1' (", dropoutRate1, ") is out of bounds [0; 1)"
)
}
if (dropoutRate2 < 0 || dropoutRate2 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate2' (", dropoutRate2, ") is out of bounds [0; 1)"
)
}
if (!is.na(eventTime) && eventTime <= 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'eventTime' (", eventTime, ") must be > 0")
}
.assertIsValidAccrualTime(accrualSetup$accrualTime)
.assertIsValidFollowUpTime(followUpTime)
.setValueAndParameterType(designPlan, "followUpTime", followUpTime, C_FOLLOW_UP_TIME_DEFAULT)
if (.isUserDefinedMaxNumberOfSubjects(designPlan) && !is.null(followUpTime) &&
length(followUpTime) == 1 && !is.na(followUpTime)) {
warning("Follow-up time will be calculated, value entered (",
followUpTime, ") is not taken into account",
call. = FALSE
)
} else if (is.na(followUpTime)) {
designPlan$followUpTime <- C_FOLLOW_UP_TIME_DEFAULT
designPlan$.setParameterType("followUpTime", C_PARAM_DEFAULT_VALUE)
}
if (objectType == "power") {
designPlan$followUpTime <- NA_real_
designPlan$.setParameterType("followUpTime", C_PARAM_NOT_APPLICABLE)
}
} else {
for (p in c(
"accrualTime", "accrualIntensity",
"eventTime", "dropoutRate1", "dropoutRate2", "dropoutTime", "followUpTime",
"analysisTime", "studyDuration"
)) {
designPlan$.setParameterType(p, C_PARAM_NOT_APPLICABLE)
}
if (designPlan$.getParameterType("accrualTime") == C_PARAM_USER_DEFINED ||
!identical(accrualTime, C_ACCRUAL_TIME_DEFAULT)) {
designPlan$.warnInCaseArgumentExists(accrualSetup$accrualTime, "accrualTime")
}
designPlan$.warnInCaseArgumentExists(dropoutRate1, "dropoutRate1")
designPlan$.warnInCaseArgumentExists(dropoutRate2, "dropoutRate2")
if (!identical(dropoutTime, C_DROP_OUT_TIME_DEFAULT)) {
designPlan$.warnInCaseArgumentExists(dropoutTime, "dropoutTime")
}
designPlan$.warnInCaseArgumentExists(maxNumberOfSubjects, "maxNumberOfSubjects")
if (!identical(followUpTime, C_FOLLOW_UP_TIME_DEFAULT)) {
designPlan$.warnInCaseArgumentExists(followUpTime, "followUpTime")
}
}
.setValueAndParameterType(designPlan, "directionUpper", directionUpper, TRUE)
if (objectType == "power") {
.setValueAndParameterType(designPlan, "maxNumberOfEvents", maxNumberOfEvents, NA_real_)
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_NOT_APPLICABLE)
}
return(designPlan)
}
.initDesignPlanSurvivalByPiecewiseSurvivalTimeObject <- function(designPlan, pwstSetup) {
designPlan$pi1 <- pwstSetup$pi1
designPlan$.setParameterType("pi1", pwstSetup$.getParameterType("pi1"))
designPlan$pi2 <- pwstSetup$pi2
designPlan$.setParameterType("pi2", pwstSetup$.getParameterType("pi2"))
designPlan$hazardRatio <- pwstSetup$hazardRatio
designPlan$.setParameterType("hazardRatio", pwstSetup$.getParameterType("hazardRatio"))
designPlan$lambda1 <- pwstSetup$lambda1
designPlan$.setParameterType("lambda1", pwstSetup$.getParameterType("lambda1"))
designPlan$lambda2 <- pwstSetup$lambda2
designPlan$.setParameterType("lambda2", pwstSetup$.getParameterType("lambda2"))
designPlan$median1 <- pwstSetup$median1
designPlan$.setParameterType("median1", pwstSetup$.getParameterType("median1"))
designPlan$median2 <- pwstSetup$median2
designPlan$.setParameterType("median2", pwstSetup$.getParameterType("median2"))
designPlan$piecewiseSurvivalTime <- pwstSetup$piecewiseSurvivalTime
designPlan$.setParameterType(
"piecewiseSurvivalTime",
pwstSetup$.getParameterType("piecewiseSurvivalTime")
)
designPlan$eventTime <- pwstSetup$eventTime
designPlan$.setParameterType("eventTime", pwstSetup$.getParameterType("eventTime"))
if (pwstSetup$.isLambdaBased()) {
return(length(designPlan$hazardRatio))
}
return(length(designPlan$pi1))
}
.initDesignPlanSurvival <- function(designPlan) {
numberOfResults <- .initDesignPlanSurvivalByPiecewiseSurvivalTimeObject(
designPlan, designPlan$.piecewiseSurvivalTime
)
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
if (length(designPlan$accountForObservationTimes) == 0 ||
is.na(designPlan$accountForObservationTimes) ||
!designPlan$accountForObservationTimes) {
designPlan$accountForObservationTimes <- TRUE
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_DEFAULT_VALUE)
}
if (!designPlan$accountForObservationTimes) {
designPlan$accountForObservationTimes <- TRUE
warning("'accountForObservationTimes' was set to TRUE ",
"because piecewise exponential survival function is enabled",
call. = FALSE
)
}
} else {
if (.isUserDefinedMaxNumberOfSubjects(designPlan)) {
if (length(designPlan$accountForObservationTimes) != 0 &&
!is.na(designPlan$accountForObservationTimes) &&
!designPlan$accountForObservationTimes) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'accountForObservationTimes' must be TRUE because 'maxNumberOfSubjects' is > 0"
)
}
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_GENERATED)
designPlan$accountForObservationTimes <- TRUE
} else {
if (length(designPlan$accountForObservationTimes) == 0 ||
is.na(designPlan$accountForObservationTimes)) {
designPlan$accountForObservationTimes <- FALSE
designPlan$.setParameterType("accountForObservationTimes", C_PARAM_DEFAULT_VALUE)
} else {
designPlan$.setParameterType(
"accountForObservationTimes",
ifelse(designPlan$accountForObservationTimes,
C_PARAM_DEFAULT_VALUE, C_PARAM_USER_DEFINED
)
)
}
}
}
designPlan$.setParameterType("chi", C_PARAM_NOT_APPLICABLE)
if (designPlan$.isSampleSizeObject()) {
designPlan$.setParameterType("directionUpper", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("maxNumberOfEvents", C_PARAM_NOT_APPLICABLE)
}
return(numberOfResults = numberOfResults)
}
.getSampleSizeFixedSurvival <- function(designPlan) {
alpha <- designPlan$getAlpha()
beta <- designPlan$getBeta()
sided <- designPlan$getSided()
twoSidedPower <- designPlan$getTwoSidedPower()
typeOfComputation <- designPlan$typeOfComputation
thetaH0 <- designPlan$thetaH0
pi1 <- designPlan$pi1
pi2 <- designPlan$pi2
allocationRatioPlanned <- designPlan$allocationRatioPlanned
accountForObservationTimes <- designPlan$accountForObservationTimes
accrualTime <- designPlan$accrualTime
kappa <- designPlan$kappa
piecewiseSurvivalTime <- designPlan$piecewiseSurvivalTime
maxNumberOfSubjects <- designPlan$maxNumberOfSubjects
hazardRatio <- designPlan$hazardRatio
.assertIsValidHazardRatio(hazardRatio, thetaH0)
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
numberOfResults <- length(hazardRatio)
} else {
numberOfResults <- length(pi1)
}
designPlan$eventsFixed <- rep(NA_real_, numberOfResults) # number of events
designPlan$nFixed <- rep(NA_real_, numberOfResults) # number of subjects
designPlan$chi <- rep(NA_real_, numberOfResults) # probability of an event
calculateAllocationRatioPlanned <- FALSE
if (allocationRatioPlanned == 0) {
allocationRatioPlannedVec <- rep(NA_real_, numberOfResults)
calculateAllocationRatioPlanned <- TRUE
designPlan$optimumAllocationRatio <- TRUE
designPlan$.setParameterType("optimumAllocationRatio", C_PARAM_USER_DEFINED)
}
userDefinedMaxNumberOfSubjects <- .isUserDefinedMaxNumberOfSubjects(designPlan)
if (userDefinedMaxNumberOfSubjects && allocationRatioPlanned == 0) {
stop(
C_EXCEPTION_TYPE_CONFLICTING_ARGUMENTS,
"determination of optimum allocation ('allocationRatioPlanned' = 0) not possible ",
"for given 'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ")"
)
}
if (userDefinedMaxNumberOfSubjects) {
timeVector <- rep(NA_real_, numberOfResults)
}
designPlan$.calculateFollowUpTime <- FALSE
lambda1 <- designPlan$lambda1
if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
lambda1 <- rep(NA_real_, numberOfResults)
}
for (i in 1:numberOfResults) {
phi <- -c(
log(1 - designPlan$dropoutRate1),
log(1 - designPlan$dropoutRate2)
) / designPlan$dropoutTime
if (!userDefinedMaxNumberOfSubjects) {
if (calculateAllocationRatioPlanned) {
# allocationRatioPlanned = 0 provides optimum sample size
allocationRatioPlanned <- stats::optimize(function(x) {
numberEvents <- .getEventsFixed(
typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
thetaH0 = thetaH0, allocationRatioPlanned = x
)
if (!accountForObservationTimes) {
probEvent <- (x * pi1[i] + pi2) / (1 + x)
} else {
probEvent <- .getEventProbabilities(
time = accrualTime[length(accrualTime)] + designPlan$followUpTime,
accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2, lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime,
phi = phi, kappa = kappa, allocationRatioPlanned = x,
hazardRatio = hazardRatio[i]
)
}
return(numberEvents / probEvent)
}, interval = c(0, 5), tol = 0.0001)$minimum
allocationRatioPlannedVec[i] <- allocationRatioPlanned
}
designPlan$eventsFixed[i] <- .getEventsFixed(
typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
thetaH0 = thetaH0, allocationRatioPlanned = allocationRatioPlanned
)
if (!accountForObservationTimes) {
designPlan$chi[i] <- (allocationRatioPlanned * pi1[i] + pi2) /
(1 + allocationRatioPlanned)
} else {
designPlan$chi[i] <- .getEventProbabilities(
time = accrualTime[length(accrualTime)] + designPlan$followUpTime,
accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio[i]
)
}
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
designPlan$nFixed[i] <- designPlan$eventsFixed[i] / designPlan$chi[i]
} else {
if (length(maxNumberOfSubjects) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of user defined 'maxNumberOfSubjects' (",
.arrayToString(maxNumberOfSubjects), ") must be 1"
)
}
designPlan$.calculateFollowUpTime <- TRUE
designPlan$eventsFixed[i] <- .getEventsFixed(
typeOfComputation = typeOfComputation, twoSidedPower = twoSidedPower,
alpha = alpha, beta = beta, sided = sided, hazardRatio = hazardRatio[i],
thetaH0 = thetaH0, allocationRatioPlanned = allocationRatioPlanned
)
designPlan$nFixed[i] <- maxNumberOfSubjects
if (designPlan$eventsFixed[i] > maxNumberOfSubjects) {
if (length(hazardRatio) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
sprintf(
paste0(
"'maxNumberOfSubjects' (%s) is smaller than the number ",
"of events (%.3f) at index %s (hazard ratio = %.3f)"
),
maxNumberOfSubjects, designPlan$eventsFixed[i], i, hazardRatio[i]
)
)
} else {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
sprintf(
paste0(
"'maxNumberOfSubjects' (%s) is smaller than the number ",
"of events (%.3f)"
),
maxNumberOfSubjects, designPlan$eventsFixed[i]
)
)
}
}
up <- 2
iterate <- 1
while (designPlan$eventsFixed[i] / .getEventProbabilities(
time = up, accrualTimeVector = accrualTime, accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2, lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = hazardRatio[i]
) > maxNumberOfSubjects) {
up <- 2 * up
iterate <- iterate + 1
if (iterate > 50) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the number of subjects is too small to reach maximum number of events ",
"(presumably due to drop-out rates), search algorithm failed"
)
}
}
timeVector[i] <- .getOneDimensionalRoot(
function(x) {
designPlan$eventsFixed[i] / .getEventProbabilities(
time = x, accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity, lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = hazardRatio[i]
) - maxNumberOfSubjects
},
lower = 0, upper = up, tolerance = 1e-06, acceptResultsOutOfTolerance = TRUE,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
if (!is.na(timeVector[i])) {
designPlan$chi[i] <- .getEventProbabilities(
time = timeVector[i],
accrualTimeVector = accrualTime,
accrualIntensity = designPlan$accrualIntensity, lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi, kappa = kappa,
allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio[i]
)
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
}
}
}
if (calculateAllocationRatioPlanned) {
allocationRatioPlanned <- allocationRatioPlannedVec
designPlan$allocationRatioPlanned <- allocationRatioPlanned
designPlan$.setParameterType("allocationRatioPlanned", C_PARAM_GENERATED)
}
if (userDefinedMaxNumberOfSubjects) {
designPlan$followUpTime <- timeVector - accrualTime[length(accrualTime)]
designPlan$.setParameterType("followUpTime", C_PARAM_GENERATED)
}
designPlan$nFixed2 <- designPlan$nFixed / (1 + allocationRatioPlanned)
designPlan$nFixed1 <- designPlan$nFixed2 * allocationRatioPlanned
designPlan$numberOfSubjects1 <- matrix(designPlan$nFixed1, nrow = 1)
designPlan$numberOfSubjects2 <- matrix(designPlan$nFixed2, nrow = 1)
if (!designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
eventRatio <- allocationRatioPlanned * pi1 / pi2
} else {
eventRatio <- NA_real_
}
# Fixed
designPlan$hazardRatio <- hazardRatio
designPlan$expectedEventsH1 <- designPlan$eventsFixed
designPlan$maxNumberOfSubjects <- designPlan$nFixed
designPlan$numberOfSubjects <- matrix(designPlan$nFixed, nrow = 1)
designPlan$.setParameterType("eventsFixed", C_PARAM_GENERATED)
designPlan$.setParameterType("nFixed1", C_PARAM_GENERATED)
designPlan$.setParameterType("nFixed2", C_PARAM_GENERATED)
designPlan$.setParameterType("nFixed", C_PARAM_GENERATED)
if (designPlan$accountForObservationTimes) {
designPlan$analysisTime <- matrix(accrualTime[length(accrualTime)] +
designPlan$followUpTime, nrow = 1)
designPlan$.setParameterType("analysisTime", C_PARAM_GENERATED)
}
return(designPlan)
}
# note that fixed sample size must be calculated before on 'designPlan'
.getSampleSizeSequentialSurvival <- function(designPlan, designCharacteristics) {
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
numberOfResults <- length(designPlan$hazardRatio)
} else {
numberOfResults <- length(designPlan$pi1)
}
kMax <- designCharacteristics$.design$kMax
designPlan$cumulativeEventsPerStage <- matrix(NA_real_, kMax, numberOfResults)
analysisTime <- matrix(NA_real_, kMax, numberOfResults)
numberOfSubjects <- matrix(NA_real_, kMax, numberOfResults)
designPlan$expectedEventsH0 <- rep(NA_real_, numberOfResults)
designPlan$expectedEventsH01 <- rep(NA_real_, numberOfResults)
designPlan$expectedEventsH1 <- rep(NA_real_, numberOfResults)
expectedNumberOfSubjectsH1 <- rep(NA_real_, numberOfResults)
studyDuration <- rep(NA_real_, numberOfResults)
designPlan$chi <- rep(NA_real_, numberOfResults)
informationRates <- designCharacteristics$information / designCharacteristics$shift
lambda1 <- designPlan$lambda1
if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
lambda1 <- rep(NA_real_, numberOfResults)
}
if (designPlan$accountForObservationTimes && designPlan$.calculateFollowUpTime) {
designPlan$followUpTime <- rep(NA_real_, numberOfResults)
}
for (i in 1:numberOfResults) {
designPlan$cumulativeEventsPerStage[, i] <- designPlan$eventsFixed[i] * informationRates *
designCharacteristics$inflationFactor
if (!designPlan$accountForObservationTimes) {
if (length(designPlan$allocationRatioPlanned) > 1) {
allocationRatioPlanned <- designPlan$allocationRatioPlanned[i]
} else {
allocationRatioPlanned <- designPlan$allocationRatioPlanned
}
designPlan$chi[i] <- (allocationRatioPlanned * designPlan$pi1[i] + designPlan$pi2) /
(1 + allocationRatioPlanned)
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
numberOfSubjects[kMax, i] <- designPlan$cumulativeEventsPerStage[kMax, i] / designPlan$chi[i]
} else {
phi <- -c(log(1 - designPlan$dropoutRate1), log(1 - designPlan$dropoutRate2)) /
designPlan$dropoutTime
if (designPlan$.calculateFollowUpTime) {
if (designPlan$cumulativeEventsPerStage[kMax, i] > designPlan$maxNumberOfSubjects[i]) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
sprintf(
paste0(
"the number of subjects (%s) is smaller than the number ",
"of events (%s) at stage %s"
),
designPlan$maxNumberOfSubjects[i],
designPlan$cumulativeEventsPerStage[kMax, i], i
)
)
}
up <- 2
iterate <- 1
while (designPlan$cumulativeEventsPerStage[kMax, i] / .getEventProbabilities(
time = up,
accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = designPlan$allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
) > designPlan$maxNumberOfSubjects[i]) {
up <- 2 * up
iterate <- iterate + 1
if (iterate > 50) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the number of subjects is too small to reach maximum number of events ",
"(presumably due to drop-out rates)"
)
}
}
totalTime <- .getOneDimensionalRoot(
function(x) {
designPlan$cumulativeEventsPerStage[kMax, i] / designPlan$maxNumberOfSubjects[i] -
.getEventProbabilities(
time = x, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = designPlan$allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = up, tolerance = 1e-06,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
# analysis times
for (j in 1:kMax) {
analysisTime[j, i] <- .getOneDimensionalRoot(
function(x) {
designPlan$cumulativeEventsPerStage[j, i] / designPlan$maxNumberOfSubjects[i] -
.getEventProbabilities(
time = x, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = designPlan$allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = totalTime, tolerance = 1e-06, acceptResultsOutOfTolerance = TRUE,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
}
analysisTime[kMax, i] <- totalTime
designPlan$followUpTime[i] <- totalTime -
designPlan$accrualTime[length(designPlan$accrualTime)]
numberOfSubjects[, i] <- .getNumberOfSubjects(
time = analysisTime[, i],
accrualTime = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = designPlan$maxNumberOfSubjects[i]
)
} else {
if (length(designPlan$allocationRatioPlanned) > 1) {
allocationRatioPlanned <- designPlan$allocationRatioPlanned[i]
} else {
allocationRatioPlanned <- designPlan$allocationRatioPlanned
}
if (is.na(designPlan$followUpTime)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'followUpTime' must be defined because 'designPlan$.calculateFollowUpTime' = FALSE"
)
}
designPlan$chi[i] <- .getEventProbabilities(
time = designPlan$accrualTime[length(designPlan$accrualTime)] + designPlan$followUpTime,
accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
designPlan$.setParameterType("chi", C_PARAM_GENERATED)
numberOfSubjects[kMax, i] <- designPlan$cumulativeEventsPerStage[kMax, i] / designPlan$chi[i]
# Analysis times
for (j in 1:(kMax - 1)) {
analysisTime[j, i] <- .getOneDimensionalRoot(
function(x) {
designPlan$cumulativeEventsPerStage[j, i] / numberOfSubjects[kMax, i] -
.getEventProbabilities(
time = x, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i],
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
phi = phi, kappa = designPlan$kappa,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = designPlan$accrualTime[length(designPlan$accrualTime)] +
designPlan$followUpTime, tolerance = 1e-06,
callingFunctionInformation = ".getSampleSizeSequentialSurvival"
)
}
analysisTime[kMax, i] <- designPlan$accrualTime[length(designPlan$accrualTime)] +
designPlan$followUpTime
numberOfSubjects[, i] <- .getNumberOfSubjects(
time = analysisTime[, i],
accrualTime = designPlan$accrualTime, accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = numberOfSubjects[kMax, i]
)
}
stoppingProbs <- designCharacteristics$rejectionProbabilities +
c(designCharacteristics$futilityProbabilities, 0)
if (all(is.na(designCharacteristics$futilityProbabilities))) {
warning("Expected number of subjects H1 and study duration H1 ",
"cannot be calculated because the futility probabilities ",
"are not applicable for the specified design",
call. = FALSE
)
}
stoppingProbs[kMax] <- 1 - sum(stoppingProbs[1:(kMax - 1)])
studyDuration[i] <- analysisTime[, i] %*% stoppingProbs
expectedNumberOfSubjectsH1[i] <- numberOfSubjects[, i] %*% stoppingProbs
}
designPlan$expectedEventsH0[i] <- designCharacteristics$averageSampleNumber0 *
designPlan$eventsFixed[i]
designPlan$expectedEventsH01[i] <- designCharacteristics$averageSampleNumber01 *
designPlan$eventsFixed[i]
designPlan$expectedEventsH1[i] <- designCharacteristics$averageSampleNumber1 *
designPlan$eventsFixed[i]
designPlan$.setParameterType("expectedEventsH0", C_PARAM_GENERATED)
designPlan$.setParameterType("expectedEventsH01", C_PARAM_GENERATED)
designPlan$.setParameterType("expectedEventsH1", C_PARAM_GENERATED)
designPlan$numberOfSubjects2 <- numberOfSubjects / (1 + designPlan$allocationRatioPlanned)
designPlan$numberOfSubjects1 <- designPlan$numberOfSubjects2 *
designPlan$allocationRatioPlanned
designPlan$.setParameterType("numberOfSubjects1", C_PARAM_GENERATED)
designPlan$.setParameterType("numberOfSubjects2", C_PARAM_GENERATED)
}
if (!is.null(designCharacteristics$rejectionProbabilities)) {
designPlan$rejectPerStage <- matrix(designCharacteristics$rejectionProbabilities,
nrow = designPlan$.design$kMax
)
designPlan$.setParameterType("rejectPerStage", C_PARAM_GENERATED)
designPlan$earlyStop <- sum(designPlan$rejectPerStage[1:(designPlan$.design$kMax - 1), ])
designPlan$.setParameterType("earlyStop", C_PARAM_GENERATED)
}
if (!is.null(designCharacteristics$futilityProbabilities) &&
any(designPlan$.design$futilityBounds != C_FUTILITY_BOUNDS_DEFAULT)) {
designPlan$futilityPerStage <- matrix(designCharacteristics$futilityProbabilities,
nrow = designPlan$.design$kMax - 1
)
designPlan$.setParameterType("futilityPerStage", C_PARAM_GENERATED)
designPlan$futilityStop <- sum(designPlan$futilityPerStage)
designPlan$.setParameterType("futilityStop", C_PARAM_GENERATED)
designPlan$earlyStop <- designPlan$earlyStop + sum(designPlan$futilityPerStage)
}
designPlan$informationRates <- matrix(informationRates, ncol = 1)
if (!is.matrix(numberOfSubjects)) {
designPlan$numberOfSubjects <- matrix(numberOfSubjects[kMax, ], nrow = 1)
} else {
designPlan$numberOfSubjects <- numberOfSubjects
}
designPlan$maxNumberOfSubjects <- designPlan$numberOfSubjects[kMax, ]
if (designPlan$.getParameterType("maxNumberOfSubjects") == C_PARAM_NOT_APPLICABLE ||
length(designPlan$maxNumberOfSubjects) > 1) {
designPlan$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
}
designPlan$maxNumberOfSubjects1 <- .getNumberOfSubjects1(
designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
)
designPlan$maxNumberOfSubjects2 <- .getNumberOfSubjects2(
designPlan$maxNumberOfSubjects, designPlan$allocationRatioPlanned
)
designPlan$.setParameterType("maxNumberOfSubjects1", C_PARAM_GENERATED)
designPlan$.setParameterType("maxNumberOfSubjects2", C_PARAM_GENERATED)
if (ncol(designPlan$informationRates) == 1 &&
identical(designPlan$informationRates[, 1], designPlan$.design$informationRates)) {
designPlan$.setParameterType("informationRates", C_PARAM_NOT_APPLICABLE)
} else {
designPlan$.setParameterType("informationRates", C_PARAM_GENERATED)
}
designPlan$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)
designPlan$.setParameterType("cumulativeEventsPerStage", C_PARAM_GENERATED)
if (designPlan$accountForObservationTimes) {
designPlan$analysisTime <- analysisTime
designPlan$expectedNumberOfSubjectsH1 <- expectedNumberOfSubjectsH1
designPlan$studyDuration <- studyDuration
designPlan$studyDurationH1 <- studyDuration # deprecated
designPlan$.setParameterType("analysisTime", C_PARAM_GENERATED)
designPlan$.setParameterType("expectedNumberOfSubjectsH1", C_PARAM_GENERATED)
designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
}
designPlan$.setParameterType("eventsFixed", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("nFixed1", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("nFixed2", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("nFixed", C_PARAM_NOT_APPLICABLE)
if (all(designPlan$allocationRatioPlanned == 1)) {
designPlan$.setParameterType("numberOfSubjects1", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("numberOfSubjects2", C_PARAM_NOT_APPLICABLE)
}
designPlan$.calculateFollowUpTime <- NA
return(designPlan)
}
.getNumberOfSubjects <- function(..., time, accrualTime, accrualIntensity, maxNumberOfSubjects) {
subjectNumbers <- c()
for (timeValue in time) {
if (is.na(timeValue)) {
return(NA_real_)
}
subjectNumbers <- c(
subjectNumbers,
.getNumberOfSubjectsInner(
timeValue = timeValue, accrualTime = accrualTime,
accrualIntensity = accrualIntensity, maxNumberOfSubjects = maxNumberOfSubjects
)
)
}
return(subjectNumbers)
}
#' @title
#' Get Event Probabilities
#'
#' @description
#' Returns the event probabilities for specified parameters at given time vector.
#'
#' @param time A numeric vector with time values.
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#' the end of accrual at specified \code{accrualIntensity} for the specified
#' number of subjects is determined or \code{accrualIntensity} is calculated
#' at fixed end of accrual.
#' @inheritParams param_three_dots
#'
#' @details
#' The function computes the overall event probabilities in a two treatment groups design.
#' For details of the parameters see \code{\link[=getSampleSizeSurvival]{getSampleSizeSurvival()}}.
#'
#' @return Returns a \code{\link{EventProbabilities}} 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.EventProbabilities]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.ParameterSet]{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
#'
#' @template examples_get_event_probabilities
#'
#' @export
#'
getEventProbabilities <- function(time, ...,
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
kappa = 1,
piecewiseSurvivalTime = NA_real_,
lambda2 = NA_real_,
lambda1 = NA_real_,
allocationRatioPlanned = 1,
hazardRatio = NA_real_,
dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
dropoutTime = 12, # C_DROP_OUT_TIME_DEFAULT
maxNumberOfSubjects = NA_real_) {
.warnInCaseOfUnknownArguments(functionName = "getEventProbabilities", ...)
.assertIsNumericVector(time, "time")
.assertIsValidMaxNumberOfSubjects(maxNumberOfSubjects, naAllowed = TRUE)
.assertIsValidAllocationRatioPlannedSampleSize(allocationRatioPlanned, maxNumberOfSubjects)
.assertIsValidKappa(kappa)
.assertIsSingleNumber(hazardRatio, "hazardRatio", naAllowed = TRUE)
if (!is.na(dropoutTime) && dropoutTime <= 0) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "'dropoutTime' (", dropoutTime, ") must be > 0", call. = FALSE)
}
if (dropoutRate1 < 0 || dropoutRate1 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate1' (", dropoutRate1, ") is out of bounds [0; 1)"
)
}
if (dropoutRate2 < 0 || dropoutRate2 >= 1) {
stop(
C_EXCEPTION_TYPE_ARGUMENT_OUT_OF_BOUNDS,
"'dropoutRate2' (", dropoutRate2, ") is out of bounds [0; 1)"
)
}
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualTime <- accrualSetup$.getAccrualTimeWithoutLeadingZero()
accrualIntensity <- accrualSetup$accrualIntensity
maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects
setting <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
hazardRatio = hazardRatio, kappa = kappa,
delayedResponseAllowed = TRUE,
.lambdaBased = TRUE
)
if (!setting$delayedResponseEnabled && length(setting$lambda1) > 1 &&
setting$.getParameterType("lambda1") == C_PARAM_USER_DEFINED) {
warning("Only the first 'lambda1' (", lambda1[1], ") was used to calculate event probabilities", call. = FALSE)
setting <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1[1],
hazardRatio = hazardRatio, kappa = kappa,
delayedResponseAllowed = TRUE,
.lambdaBased = TRUE
)
}
piecewiseSurvivalTime <- setting$piecewiseSurvivalTime
lambda2 <- setting$lambda2
lambda1 <- setting$lambda1
hazardRatio <- setting$hazardRatio
phi <- -log(1 - c(dropoutRate1, dropoutRate2)) / dropoutTime
if (length(accrualTime) != length(accrualIntensity)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of 'accrualTime' (", (length(accrualTime) + 1),
") must be equal to length of 'accrualIntensity' (", length(accrualIntensity), ")"
)
}
if (any(accrualIntensity <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualIntensity' must be > 0")
}
if (any(accrualTime <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualTime' must be > 0")
}
if (kappa != 1 && any(phi > 0)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"for Weibull distribution (kappa != 1) drop-out rates (phi) cannot be specified"
)
}
if (any(phi < 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all drop-out rates (phi) must be >= 0")
}
.assertIsNumericVector(lambda2, "lambda2")
if (any(lambda2 <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all rates (lambda2) must be > 0")
}
eventProbabilities <- EventProbabilities(
.piecewiseSurvivalTime = setting,
.accrualTime = accrualSetup,
time = time,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda1 = lambda1,
lambda2 = lambda2,
allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = hazardRatio,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
maxNumberOfSubjects = maxNumberOfSubjects
)
eventProbabilities$.setParameterType("time", C_PARAM_USER_DEFINED)
eventProbabilities$.setParameterType(
"accrualTime",
accrualSetup$.getParameterType("accrualTime")
)
eventProbabilities$.setParameterType(
"accrualIntensity",
accrualSetup$.getParameterType("accrualIntensity")
)
eventProbabilities$.setParameterType("kappa", setting$.getParameterType("kappa"))
eventProbabilities$.setParameterType(
"piecewiseSurvivalTime",
setting$.getParameterType("piecewiseSurvivalTime")
)
eventProbabilities$.setParameterType("lambda1", setting$.getParameterType("lambda1"))
eventProbabilities$.setParameterType("lambda2", setting$.getParameterType("lambda2"))
.setValueAndParameterType(eventProbabilities, "allocationRatioPlanned", allocationRatioPlanned, 1)
eventProbabilities$.setParameterType("hazardRatio", setting$.getParameterType("hazardRatio"))
.setValueAndParameterType(eventProbabilities, "dropoutRate1", dropoutRate1, C_DROP_OUT_RATE_1_DEFAULT)
.setValueAndParameterType(eventProbabilities, "dropoutRate2", dropoutRate2, C_DROP_OUT_RATE_2_DEFAULT)
.setValueAndParameterType(eventProbabilities, "dropoutTime", dropoutTime, C_DROP_OUT_TIME_DEFAULT)
eventProbabilities$.setParameterType(
"maxNumberOfSubjects",
accrualSetup$.getParameterType("maxNumberOfSubjects")
)
eventProbabilities$cumulativeEventProbabilities <- numeric(0)
eventProbabilities$overallEventProbabilities <- numeric(0) # deprecated
eventProbabilities$eventProbabilities1 <- numeric(0)
eventProbabilities$eventProbabilities2 <- numeric(0)
for (timeValue in time) {
eventProbs <- .getEventProbabilitiesGroupwise(
time = timeValue,
accrualTimeVector = accrualSetup$.getAccrualTimeWithoutLeadingZero(),
accrualIntensity = accrualSetup$accrualIntensity, lambda2 = lambda2,
lambda1 = lambda1, piecewiseSurvivalTime = piecewiseSurvivalTime, phi = phi,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned, hazardRatio = hazardRatio
)
eventProbabilities$cumulativeEventProbabilities <- c(
eventProbabilities$cumulativeEventProbabilities,
.getEventProbabilitiesOverall(eventProbs, allocationRatioPlanned)
)
eventProbabilities$overallEventProbabilities <-
eventProbabilities$cumulativeEventProbabilities # deprecated
eventProbabilities$eventProbabilities1 <- c(
eventProbabilities$eventProbabilities1,
eventProbs[1]
)
eventProbabilities$eventProbabilities2 <- c(
eventProbabilities$eventProbabilities2,
eventProbs[2]
)
}
eventProbabilities$.setParameterType("cumulativeEventProbabilities", C_PARAM_GENERATED)
eventProbabilities$.setParameterType("eventProbabilities1", C_PARAM_GENERATED)
eventProbabilities$.setParameterType("eventProbabilities2", C_PARAM_GENERATED)
return(eventProbabilities)
}
#' @title
#' Get Number Of Subjects
#'
#' @description
#' Returns the number of recruited subjects at given time vector.
#'
#' @param time A numeric vector with time values.
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#' the end of accrual at specified \code{accrualIntensity} for the specified number of
#' subjects is determined or \code{accrualIntensity} is calculated at fixed end of accrual.
#' @inheritParams param_three_dots
#'
#' @details
#' Calculate number of subjects over time range at given accrual time vector
#' and accrual intensity. Intensity can either be defined in absolute or
#' relative terms (for the latter, \code{maxNumberOfSubjects} needs to be defined)\cr
#' The function is used by \code{\link[=getSampleSizeSurvival]{getSampleSizeSurvival()}}.
#'
#' @return Returns a \code{\link{NumberOfSubjects}} 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.NumberOfSubjects]{plot()}} to plot the object,
#' \item \code{\link[=as.data.frame.ParameterSet]{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
#'
#' @seealso \code{\link{AccrualTime}} for defining the accrual time.
#'
#' @template examples_get_number_of_subjects
#'
#' @export
#'
getNumberOfSubjects <- function(time, ...,
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
maxNumberOfSubjects = NA_real_) {
.warnInCaseOfUnknownArguments(functionName = "getNumberOfSubjects", ...)
.assertIsNumericVector(time, "time")
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects
)
accrualTime <- accrualSetup$.getAccrualTimeWithoutLeadingZero()
accrualIntensity <- accrualSetup$accrualIntensity
maxNumberOfSubjects <- accrualSetup$maxNumberOfSubjects
if (length(accrualTime) != length(accrualIntensity)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"length of 'accrualTime' (", length(accrualTime),
") must be equal to length of 'accrualIntensity' (", length(accrualIntensity), ")"
)
}
if (any(accrualIntensity < 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualIntensity' must be >= 0")
}
if (all(accrualIntensity < 1)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "at least one value of 'accrualIntensity' must be >= 1")
}
if (any(accrualTime <= 0)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT, "all values of 'accrualTime' must be > 0")
}
numberOfSubjects <- .getNumberOfSubjects(
time = time, accrualTime = accrualTime,
accrualIntensity = accrualIntensity, maxNumberOfSubjects = maxNumberOfSubjects
)
result <- NumberOfSubjects(
.accrualTime = accrualSetup,
time = time,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
maxNumberOfSubjects = maxNumberOfSubjects,
numberOfSubjects = numberOfSubjects
)
result$.setParameterType("time", C_PARAM_USER_DEFINED)
result$.setParameterType("accrualTime", accrualSetup$.getParameterType("accrualTime"))
result$.setParameterType("accrualIntensity", accrualSetup$.getParameterType("accrualIntensity"))
result$.setParameterType("maxNumberOfSubjects", accrualSetup$.getParameterType("maxNumberOfSubjects"))
result$.setParameterType("numberOfSubjects", C_PARAM_GENERATED)
return(result)
}
#' @title
#' Get Sample Size Survival
#'
#' @description
#' Returns the sample size for testing the hazard ratio in a two treatment groups survival design.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_typeOfComputation
#' @inheritParams param_allocationRatioPlanned_sampleSize
#' @inheritParams param_thetaH0
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_pi1_survival
#' @inheritParams param_pi2_survival
#' @inheritParams param_median1
#' @inheritParams param_median2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_eventTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param followUpTime The assumed (additional) follow-up time for the study, default is \code{6}.
#' The total study duration is \code{accrualTime + followUpTime}.
#' @param maxNumberOfSubjects If \code{maxNumberOfSubjects > 0} is specified,
#' the follow-up time for the required number of events is determined.
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the number of events and an estimate for the
#' necessary number of subjects for testing the hazard ratio in a survival design.
#' It also calculates the time when the required events are expected under the given
#' assumptions (exponentially, piecewise exponentially, or Weibull distributed survival times
#' and constant or non-constant piecewise accrual).
#' Additionally, an allocation ratio = \code{n1 / n2} can be specified where \code{n1} and \code{n2} are the number
#' of subjects in the two treatment groups.
#'
#' Optional argument \code{accountForObservationTimes}: if \code{accountForObservationTimes = TRUE}, the number of
#' subjects is calculated assuming specific accrual and follow-up time, default is \code{TRUE}.
#'
#' The formula of Kim & Tsiatis (Biometrics, 1990)
#' is used to calculate the expected number of events under the alternative
#' (see also Lakatos & Lan, Statistics in Medicine, 1992). These formulas are generalized
#' to piecewise survival times and non-constant piecewise accrual over time.\cr
#'
#' Optional argument \code{accountForObservationTimes}: if \code{accountForObservationTimes = FALSE},
#' only the event rates are used for the calculation of the maximum number of subjects.
#'
#' @template details_piecewise_survival
#'
#' @template details_piecewise_accrual
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family sample size functions
#'
#' @template examples_get_sample_size_survival
#'
#' @export
#'
getSampleSizeSurvival <- function(design = NULL, ...,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0 = 1, # C_THETA_H0_SURVIVAL_DEFAULT
pi1 = NA_real_,
pi2 = NA_real_,
lambda1 = NA_real_,
lambda2 = NA_real_,
median1 = NA_real_,
median2 = NA_real_,
kappa = 1,
hazardRatio = NA_real_,
piecewiseSurvivalTime = NA_real_,
allocationRatioPlanned = NA_real_, # C_ALLOCATION_RATIO_DEFAULT
eventTime = 12, # C_EVENT_TIME_DEFAULT
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
followUpTime = NA_real_,
maxNumberOfSubjects = NA_real_,
dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
dropoutTime = 12 # C_DROP_OUT_TIME_DEFAULT
) {
if (is.null(design)) {
design <- .getDefaultDesign(..., type = "sampleSize", ignore = c("accountForObservationTimes"))
.warnInCaseOfUnknownArguments(
functionName = "getSampleSizeSurvival",
ignore = c(.getDesignArgumentsToIgnoreAtUnknownArgumentCheck(
design,
powerCalculationEnabled = FALSE
), "accountForObservationTimes"), ...
)
} else {
.assertIsTrialDesign(design)
.warnInCaseOfUnknownArguments(
functionName = "getSampleSizeSurvival", ...,
ignore = c("accountForObservationTimes")
)
.warnInCaseOfTwoSidedPowerArgument(...)
}
if (!is.na(maxNumberOfSubjects) && maxNumberOfSubjects == 0) {
maxNumberOfSubjects <- NA_real_
}
# identify accrual time case
accrualSetup <- getAccrualTime(
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
maxNumberOfSubjects = maxNumberOfSubjects,
showWarnings = FALSE
)
accrualSetup$.validate()
accountForObservationTimes <- .getOptionalArgument("accountForObservationTimes", ...)
if (is.null(accountForObservationTimes)) {
accountForObservationTimes <- TRUE
}
if (!accrualSetup$maxNumberOfSubjectsCanBeCalculatedDirectly &&
accrualSetup$followUpTimeMustBeUserDefined) {
if (is.na(followUpTime)) {
if (accrualSetup$piecewiseAccrualEnabled && !accrualSetup$endOfAccrualIsUserDefined) {
stop(
C_EXCEPTION_TYPE_MISSING_ARGUMENT,
"'followUpTime', 'maxNumberOfSubjects' or end of accrual must be defined"
)
}
stop(
C_EXCEPTION_TYPE_MISSING_ARGUMENT,
"'followUpTime' or 'maxNumberOfSubjects' must be defined"
)
}
if (followUpTime == Inf) {
followUpTime <- 1e12
}
if (!any(is.na(hazardRatio)) && !is.na(thetaH0)) {
.assertIsValidHazardRatio(hazardRatio, thetaH0)
}
pwst <- getPiecewiseSurvivalTime(
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda1 = lambda1, lambda2 = lambda2,
pi1 = pi1, pi2 = pi2,
median1 = median1, median2 = median2,
hazardRatio = hazardRatio, eventTime = eventTime, kappa = kappa,
.silent = TRUE
)
paramName <- NULL
if (!pwst$piecewiseSurvivalEnabled) {
if (pwst$.getParameterType("pi1") == C_PARAM_USER_DEFINED ||
pwst$.getParameterType("pi1") == C_PARAM_DEFAULT_VALUE ||
pwst$.getParameterType("pi2") == C_PARAM_USER_DEFINED) {
paramName <- "pi1"
} else if (pwst$.getParameterType("lambda1") == C_PARAM_USER_DEFINED ||
pwst$.getParameterType("lambda2") == C_PARAM_USER_DEFINED) {
paramName <- "lambda1"
} else if (pwst$.getParameterType("hazardRatio") == C_PARAM_USER_DEFINED) {
paramName <- "hazardRatio"
} else if (pwst$.getParameterType("median1") == C_PARAM_USER_DEFINED ||
pwst$.getParameterType("median2") == C_PARAM_USER_DEFINED) {
paramName <- "median1"
}
} else if (pwst$.getParameterType("hazardRatio") == C_PARAM_USER_DEFINED) {
paramName <- "hazardRatio"
}
if (!is.null(paramName)) {
paramValue <- pwst[[paramName]]
if (!is.null(paramValue) && length(paramValue) > 1) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"the calculation of 'maxNumberOfSubjects' for given 'followUpTime' ",
"is only available for a single '", paramName, "'; ",
paramName, " = ", .arrayToString(
paramValue,
vectorLookAndFeelEnabled = TRUE
)
)
}
}
hr <- hazardRatio
if (all(is.na(hazardRatio))) {
hr <- pwst$hazardRatio
}
if (all(is.na(hazardRatio))) {
.assertIsValidHazardRatio(hr, thetaH0)
}
maxNumberOfSubjectsTarget <- NA_real_
withCallingHandlers(
{
# search for accrual time that provides a result
at <- accrualSetup$accrualTime
additionalAccrual <- 1
searchAccrualTimeEnabled <- TRUE
maxSearchIterations <- 50
maxNumberOfSubjectsLower <- NA_real_
maxNumberOfSubjectsLowerBefore <- 0
sampleSize <- NULL
expectionMessage <- NA_character_
while (searchAccrualTimeEnabled && maxSearchIterations >= 0 &&
(is.na(maxNumberOfSubjectsLower) ||
maxNumberOfSubjectsLower < maxNumberOfSubjectsLowerBefore ||
maxNumberOfSubjectsLower < 1e8)) {
tryCatch(
{
maxNumberOfSubjectsLowerBefore <- ifelse(is.na(maxNumberOfSubjectsLower),
0, maxNumberOfSubjectsLower
)
maxNumberOfSubjectsLower <- getAccrualTime(
accrualTime = c(at, at[length(at)] + additionalAccrual),
accrualIntensity = accrualSetup$accrualIntensity,
accrualIntensityType = accrualIntensityType
)$maxNumberOfSubjects
additionalAccrual <- 2 * additionalAccrual
sampleSize <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1, median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsLower,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2, dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)
searchAccrualTimeEnabled <- FALSE
},
error = function(e) {
expectionMessage <<- e$message
}
)
maxSearchIterations <- maxSearchIterations - 1
}
if (is.null(sampleSize) || is.na(sampleSize$followUpTime)) {
if (!is.na(expectionMessage) && grepl("'allocationRatioPlanned' > 0", expectionMessage)) {
stop(expectionMessage, call. = FALSE)
}
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"'additionalAccrual' could not be found, change accrual time specification",
call. = FALSE
)
}
# define lower bound for maxNumberOfSubjects
maxNumberOfSubjectsLower <- ceiling(max(na.omit(c(
sampleSize$eventsFixed,
as.vector(sampleSize$cumulativeEventsPerStage)
))))
if (is.na(maxNumberOfSubjectsLower)) {
stop(C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'maxNumberOfSubjectsLower' could not be found",
call. = FALSE
)
}
# check whether accrual time already fulfills requirement
# (followUpTime < given value) or need to be increased,
# then define upper bound for maxNumberOfSubjects
maxSearchIterations <- 50
maxNumberOfSubjectsUpper <- NA_real_
fut <- sampleSize$followUpTime
iterations <- 1
while (fut <= followUpTime) {
fut <- 2 * abs(fut)
iterations <- iterations + 1
if (iterations > 50) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"search algorithm failed to end",
call. = FALSE
)
}
}
while (!is.na(fut) && fut > followUpTime && maxSearchIterations >= 0) {
maxNumberOfSubjectsUpper <- getAccrualTime(
accrualTime = c(at, at[length(at)] + additionalAccrual),
accrualIntensity = accrualSetup$accrualIntensity,
accrualIntensityType = accrualIntensityType
)$maxNumberOfSubjects
additionalAccrual <- 2 * additionalAccrual
fut <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1, median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsUpper,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2, dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)$followUpTime
maxSearchIterations <- maxSearchIterations - 1
}
if (is.na(maxNumberOfSubjectsUpper)) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"'maxNumberOfSubjectsUpper' could not be found ",
"(fut = ", fut, ", followUpTime = ", followUpTime, ")",
call. = FALSE
)
}
# use maxNumberOfSubjectsLower and maxNumberOfSubjectsUpper to find end of accrual
if (dropoutRate1 != 0 || dropoutRate2 != 0) {
# Adjust lower bound for given dropouts assuming exponential distribution
if (is.na(allocationRatioPlanned)) {
allocationRatioPlanned <- C_ALLOCATION_RATIO_DEFAULT
}
maxNumberOfSubjectsLower <- maxNumberOfSubjectsLower /
((allocationRatioPlanned * (1 - dropoutRate1)^(
accrualSetup$accrualTime[length(accrualSetup$accrualTime)] / dropoutTime) +
(1 - dropoutRate2)^(accrualSetup$accrualTime[length(accrualSetup$accrualTime)] / dropoutTime)) /
(allocationRatioPlanned + 1))
prec <- 1
maxSearchIterations <- 50
while (prec > 1e-04 && maxSearchIterations >= 0) {
maxNumberOfSubjectsTarget <- (maxNumberOfSubjectsLower + maxNumberOfSubjectsUpper) / 2
fut <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsTarget,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = FALSE
)$followUpTime
ifelse(fut <= followUpTime,
maxNumberOfSubjectsUpper <- maxNumberOfSubjectsTarget,
maxNumberOfSubjectsLower <- maxNumberOfSubjectsTarget
)
prec <- maxNumberOfSubjectsUpper - maxNumberOfSubjectsLower
maxSearchIterations <- maxSearchIterations - 1
}
} else {
maxNumberOfSubjectsTarget <- .getOneDimensionalRootBisectionMethod(
fun = function(x) {
fut <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = x,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = FALSE
)$followUpTime
return(followUpTime - fut)
},
lower = maxNumberOfSubjectsLower,
upper = maxNumberOfSubjectsUpper,
tolerance = 1e-04,
acceptResultsOutOfTolerance = TRUE,
maxSearchIterations = 50,
direction = 0,
suppressWarnings = FALSE,
callingFunctionInformation = "getSampleSizeSurvival"
)
}
},
warning = function(w) {
invokeRestart("muffleWarning")
}
)
if (is.na(maxNumberOfSubjectsTarget)) {
stop(
C_EXCEPTION_TYPE_RUNTIME_ISSUE,
"failed to calculate 'maxNumberOfSubjects' by given 'followUpTime' ",
"(lower = ", maxNumberOfSubjectsLower, ", upper = ", maxNumberOfSubjectsUpper, ")"
)
}
sampleSizeSurvival <- .getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualSetup$accrualTime,
accrualIntensity = accrualSetup$accrualIntensity, kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime, lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = NA_real_, maxNumberOfSubjects = maxNumberOfSubjectsTarget,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio,
userFunctionCallEnabled = FALSE
)
sampleSizeSurvival$.setParameterType("followUpTime", C_PARAM_USER_DEFINED)
sampleSizeSurvival$.accrualTime <- accrualSetup
if (!is.na(sampleSizeSurvival$followUpTime)) {
if (followUpTime == 1e12) {
followUpTime <- Inf
}
if (sampleSizeSurvival$followUpTime >= -1e-02 && sampleSizeSurvival$followUpTime <= 1e-02) {
sampleSizeSurvival$followUpTime <- 0
}
if (sampleSizeSurvival$followUpTime < followUpTime - 1e-02 ||
sampleSizeSurvival$followUpTime > followUpTime + 1e-02) {
sampleSizeSurvival$.setParameterType("followUpTime", C_PARAM_GENERATED)
warning("User defined 'followUpTime' (", followUpTime, ") ignored because ",
"follow-up time is ", round(sampleSizeSurvival$followUpTime, 4),
call. = FALSE
)
}
}
sampleSizeSurvival$.setParameterType("maxNumberOfSubjects", C_PARAM_GENERATED)
sampleSizeSurvival$.setParameterType("accrualTime", C_PARAM_GENERATED)
return(sampleSizeSurvival)
}
return(.getSampleSizeSurvival(
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0, pi2 = pi2, pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = accountForObservationTimes,
eventTime = eventTime, accrualTime = accrualTime,
accrualIntensity = accrualIntensity, accrualIntensityType = accrualIntensityType,
kappa = kappa, piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2, lambda1 = lambda1,
median1 = median1, median2 = median2,
followUpTime = followUpTime, maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1, dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio
))
}
#' @title
#' Get Power Survival
#'
#' @description
#' Returns the power, stopping probabilities, and expected sample size for testing
#' the hazard ratio in a two treatment groups survival design.
#'
#' @inheritParams param_design_with_default
#' @inheritParams param_typeOfComputation
#' @inheritParams param_allocationRatioPlanned
#' @inheritParams param_thetaH0
#' @inheritParams param_lambda1
#' @inheritParams param_lambda2
#' @inheritParams param_pi1_survival
#' @inheritParams param_pi2_survival
#' @inheritParams param_median1
#' @inheritParams param_median2
#' @inheritParams param_piecewiseSurvivalTime
#' @inheritParams param_directionUpper
#' @inheritParams param_accrualTime
#' @inheritParams param_accrualIntensity
#' @inheritParams param_accrualIntensityType
#' @inheritParams param_eventTime
#' @inheritParams param_hazardRatio
#' @inheritParams param_kappa
#' @inheritParams param_dropoutRate1
#' @inheritParams param_dropoutRate2
#' @inheritParams param_dropoutTime
#' @param maxNumberOfEvents \code{maxNumberOfEvents > 0} is the maximum number of events, it determines
#' the power of the test and needs to be specified.
#' @inheritParams param_maxNumberOfSubjects_survival
#' @inheritParams param_three_dots
#'
#' @details
#' At given design the function calculates the power, stopping probabilities, and expected
#' sample size at given number of events and number of subjects.
#' It also calculates the time when the required events are expected under the given
#' assumptions (exponentially, piecewise exponentially, or Weibull distributed survival times
#' and constant or non-constant piecewise accrual).
#' Additionally, an allocation ratio = \code{n1 / n2} can be specified where \code{n1} and \code{n2} are the number
#' of subjects in the two treatment groups.
#'
#' The formula of Kim & Tsiatis (Biometrics, 1990)
#' is used to calculate the expected number of events under the alternative
#' (see also Lakatos & Lan, Statistics in Medicine, 1992). These formulas are generalized to piecewise survival times and
#' non-constant piecewise accrual over time.\cr
#'
#' @template details_piecewise_survival
#'
#' @template details_piecewise_accrual
#'
#' @template return_object_trial_design_plan
#' @template how_to_get_help_for_generics
#'
#' @family power functions
#'
#' @template examples_get_power_survival
#'
#' @export
#'
getPowerSurvival <- function(design = NULL, ...,
typeOfComputation = c("Schoenfeld", "Freedman", "HsiehFreedman"),
thetaH0 = 1, # C_THETA_H0_SURVIVAL_DEFAULT
directionUpper = NA,
pi1 = NA_real_,
pi2 = NA_real_,
lambda1 = NA_real_,
lambda2 = NA_real_,
median1 = NA_real_,
median2 = NA_real_,
kappa = 1,
hazardRatio = NA_real_,
piecewiseSurvivalTime = NA_real_,
allocationRatioPlanned = 1, # C_ALLOCATION_RATIO_DEFAULT
eventTime = 12, # C_EVENT_TIME_DEFAULT
accrualTime = c(0, 12), # C_ACCRUAL_TIME_DEFAULT
accrualIntensity = 0.1, # C_ACCRUAL_INTENSITY_DEFAULT
accrualIntensityType = c("auto", "absolute", "relative"),
maxNumberOfSubjects = NA_real_,
maxNumberOfEvents = NA_real_,
dropoutRate1 = 0, # C_DROP_OUT_RATE_1_DEFAULT
dropoutRate2 = 0, # C_DROP_OUT_RATE_2_DEFAULT
dropoutTime = 12 # C_DROP_OUT_TIME_DEFAULT
) {
if (is.null(design)) {
design <- .getDefaultDesign(..., type = "power")
.warnInCaseOfUnknownArguments(
functionName = "getPowerSurvival",
ignore = .getDesignArgumentsToIgnoreAtUnknownArgumentCheck(design, powerCalculationEnabled = TRUE), ...
)
} else {
.assertIsTrialDesign(design)
.warnInCaseOfUnknownArguments(functionName = "getPowerSurvival", ...)
.warnInCaseOfTwoSidedPowerArgument(...)
.warnInCaseOfTwoSidedPowerIsDisabled(design)
}
designPlan <- .createDesignPlanSurvival(
objectType = "power",
design = design,
typeOfComputation = typeOfComputation,
thetaH0 = thetaH0,
pi2 = pi2,
pi1 = pi1,
allocationRatioPlanned = allocationRatioPlanned,
accountForObservationTimes = TRUE,
eventTime = eventTime,
accrualTime = accrualTime,
accrualIntensity = accrualIntensity,
accrualIntensityType = accrualIntensityType,
kappa = kappa,
piecewiseSurvivalTime = piecewiseSurvivalTime,
lambda2 = lambda2,
lambda1 = lambda1,
median1 = median1,
median2 = median2,
directionUpper = directionUpper,
maxNumberOfEvents = maxNumberOfEvents,
maxNumberOfSubjects = maxNumberOfSubjects,
dropoutRate1 = dropoutRate1,
dropoutRate2 = dropoutRate2,
dropoutTime = dropoutTime,
hazardRatio = hazardRatio
)
.assertIsSingleNumber(allocationRatioPlanned, "allocationRatioPlanned")
.assertIsInOpenInterval(allocationRatioPlanned, "allocationRatioPlanned", 0, C_ALLOCATION_RATIO_MAXIMUM)
if (designPlan$typeOfComputation == "Schoenfeld") {
theta <- sqrt(allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(log(designPlan$hazardRatio / thetaH0))
} else if (designPlan$typeOfComputation == "Freedman") {
theta <- sqrt(allocationRatioPlanned) * (designPlan$hazardRatio - 1) /
(allocationRatioPlanned * designPlan$hazardRatio + 1)
} else if (designPlan$typeOfComputation == "HsiehFreedman") {
theta <- sqrt(4 * allocationRatioPlanned) / (1 + allocationRatioPlanned) *
(designPlan$hazardRatio - 1) / (designPlan$hazardRatio + 1)
}
if (!is.na(designPlan$directionUpper) && !designPlan$directionUpper) {
theta <- -theta
}
powerAndAverageSampleNumber <- getPowerAndAverageSampleNumber(
design = design, theta = theta, nMax = maxNumberOfEvents
)
kMax <- design$kMax
sided <- design$sided
if (designPlan$.piecewiseSurvivalTime$.isLambdaBased()) {
numberOfResults <- length(designPlan$hazardRatio)
} else {
numberOfResults <- length(designPlan$pi1)
}
stoppingProbs <- matrix(NA_real_, kMax, numberOfResults)
designPlan$analysisTime <- matrix(NA_real_, kMax, numberOfResults)
designPlan$numberOfSubjects <- matrix(NA_real_, kMax, numberOfResults)
designPlan$studyDuration <- rep(NA_real_, numberOfResults)
designPlan$expectedNumberOfSubjects <- rep(NA_real_, numberOfResults)
cumulativeEventsPerStage <- maxNumberOfEvents * design$informationRates
parameterNames <- c(
"analysisTime",
"numberOfSubjects",
"studyDuration",
"expectedNumberOfSubjects",
"cumulativeEventsPerStage"
)
phi <- -c(log(1 - dropoutRate1), log(1 - dropoutRate2)) / dropoutTime
lambda1 <- designPlan$lambda1
if (designPlan$.piecewiseSurvivalTime$piecewiseSurvivalEnabled) {
lambda1 <- rep(NA_real_, numberOfResults)
}
for (i in 1:numberOfResults) {
# Analysis times
up <- 2
iterate <- 1
while (cumulativeEventsPerStage[kMax] / designPlan$maxNumberOfSubjects > .getEventProbabilities(
time = up, accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i], phi = phi,
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)) {
up <- 2 * up
iterate <- iterate + 1
if (iterate > 50) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ") ",
"is too small to reach maximum number of events ",
"(presumably due to drop-out rates)"
)
}
}
for (j in 1:kMax) {
designPlan$analysisTime[j, i] <- .getOneDimensionalRoot(
function(x) {
cumulativeEventsPerStage[j] / designPlan$maxNumberOfSubjects -
.getEventProbabilities(
time = x,
accrualTimeVector = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
lambda2 = designPlan$lambda2,
lambda1 = lambda1[i], phi = phi,
piecewiseSurvivalTime = designPlan$piecewiseSurvivalTime,
kappa = kappa, allocationRatioPlanned = allocationRatioPlanned,
hazardRatio = designPlan$hazardRatio[i]
)
},
lower = 0, upper = up, tolerance = 1e-06,
callingFunctionInformation = "getPowerSurvival"
)
if (is.na(designPlan$analysisTime[j, i])) {
warning("Cannot calculate analysis time at stage ", j, ": ",
"'maxNumberOfSubjects' (", designPlan$maxNumberOfSubjects, ") is too ",
"small to reach maximum number of events",
call. = FALSE
)
}
}
if (kMax > 1) {
designPlan$numberOfSubjects[, i] <- .getNumberOfSubjects(
time = designPlan$analysisTime[, i],
accrualTime = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = designPlan$maxNumberOfSubjects
)
powerAndAverageSampleNumber$futilityPerStage[is.na(
powerAndAverageSampleNumber$futilityPerStage[, i]
), i] <- 0
stoppingProbs[, i] <- powerAndAverageSampleNumber$rejectPerStage[, i] +
c(powerAndAverageSampleNumber$futilityPerStage[, i], 0)
stoppingProbs[kMax, i] <- 1 - sum(stoppingProbs[1:(kMax - 1), i])
designPlan$studyDuration[i] <- designPlan$analysisTime[, i] %*% stoppingProbs[, i]
designPlan$.setParameterType("studyDuration", C_PARAM_GENERATED)
designPlan$expectedNumberOfSubjects[i] <-
designPlan$numberOfSubjects[, i] %*% stoppingProbs[, i]
}
}
if (kMax == 1) {
designPlan$expectedNumberOfSubjects <- .getNumberOfSubjects(
time = designPlan$analysisTime[1, ],
accrualTime = designPlan$accrualTime,
accrualIntensity = designPlan$accrualIntensity,
maxNumberOfSubjects = designPlan$maxNumberOfSubjects
)
designPlan$numberOfSubjects <- matrix(designPlan$expectedNumberOfSubjects, nrow = 1)
}
designPlan$cumulativeEventsPerStage <- matrix(cumulativeEventsPerStage, ncol = 1)
designPlan$.setParameterType("cumulativeEventsPerStage", C_PARAM_GENERATED)
.addDeprecatedFieldValues(designPlan, "eventsPerStage", designPlan$cumulativeEventsPerStage)
designPlan$expectedNumberOfEvents <- powerAndAverageSampleNumber$averageSampleNumber
designPlan$overallReject <- powerAndAverageSampleNumber$overallReject
designPlan$rejectPerStage <- powerAndAverageSampleNumber$rejectPerStage
designPlan$futilityStop <- powerAndAverageSampleNumber$overallFutility
designPlan$futilityPerStage <- powerAndAverageSampleNumber$futilityPerStage
designPlan$earlyStop <- powerAndAverageSampleNumber$overallEarlyStop
parameterNames <- c(
parameterNames,
"expectedNumberOfEvents",
"overallReject",
"rejectPerStage",
"futilityStop",
"futilityPerStage",
"earlyStop"
)
for (parameterName in parameterNames) {
designPlan$.setParameterType(parameterName, C_PARAM_GENERATED)
}
if (kMax == 1L) {
designPlan$.setParameterType("numberOfSubjects", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("cumulativeEventsPerStage", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("futilityStop", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("earlyStop", C_PARAM_NOT_APPLICABLE)
designPlan$.setParameterType("rejectPerStage", C_PARAM_NOT_APPLICABLE)
}
if (!any(is.na(designPlan$analysisTime)) && !any(is.na(designPlan$accrualTime))) {
designPlan$followUpTime <- designPlan$analysisTime[kMax, ] -
designPlan$accrualTime[length(designPlan$accrualTime)]
designPlan$.setParameterType("followUpTime", C_PARAM_GENERATED)
}
designPlan$singleEventsPerStage <- .getSingleEventsPerStage(designPlan$cumulativeEventsPerStage)
.addEffectScaleBoundaryDataToDesignPlan(designPlan)
.addStudyDurationToDesignPlan(designPlan)
.hideFutilityStopsIfNotApplicable(designPlan)
return(designPlan)
}
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.