Nothing
## |
## | *Simulation of enrichment design with binary data*
## |
## | 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: 7383 $
## | Last changed: $Date: 2023-11-02 15:18:21 +0100 (Do, 02 Nov 2023) $
## | Last changed by: $Author: pahlke $
## |
#' @include f_simulation_enrichment.R
NULL
.getSimulationRatesEnrichmentStageSubjects <- function(...,
stage,
directionUpper,
conditionalPower,
conditionalCriticalValue,
plannedSubjects,
allocationRatioPlanned,
selectedPopulations,
piTreatmentH1,
piControlH1,
overallRatesTreatment,
overallRatesControl,
minNumberOfSubjectsPerStage,
maxNumberOfSubjectsPerStage) {
stage <- stage - 1 # to be consistent with non-enrichment situation
gMax <- nrow(overallRatesTreatment)
if (!is.na(conditionalPower)) {
if (any(selectedPopulations[1:gMax, stage + 1], na.rm = TRUE)) {
if (is.na(piControlH1)) {
pi2H1 <- overallRatesControl[selectedPopulations[1:gMax, stage + 1], stage]
} else {
pi2H1 <- piControlH1
}
if (is.na(piTreatmentH1)) {
pi1H1 <- overallRatesTreatment[selectedPopulations[1:gMax, stage + 1], stage]
} else {
pi1H1 <- piTreatmentH1
}
pim <- (allocationRatioPlanned[stage] * pi1H1 + pi2H1) / (1 + allocationRatioPlanned[stage])
if (conditionalCriticalValue[stage] > 8) {
newSubjects <- maxNumberOfSubjectsPerStage[stage + 1]
} else {
newSubjects <- (1 + 1 / allocationRatioPlanned[stage]) * (max(0, conditionalCriticalValue[stage] *
sqrt(pim * (1 - pim) * (1 + allocationRatioPlanned[stage])) +
.getQNorm(conditionalPower) * sqrt(pi1H1 * (1 - pi1H1) +
pi2H1 * (1 - pi2H1) * allocationRatioPlanned[stage]), na.rm = TRUE))^2 /
(max(1e-07, (2 * directionUpper - 1) * (pi1H1 - pi2H1), na.rm = TRUE))^2
newSubjects <- min(
max(minNumberOfSubjectsPerStage[stage + 1], newSubjects),
maxNumberOfSubjectsPerStage[stage + 1]
)
}
} else {
newSubjects <- 0
}
} else {
newSubjects <- plannedSubjects[stage + 1] - plannedSubjects[stage]
}
return(newSubjects)
}
.getSimulatedStageRatesEnrichment <- function(...,
design,
subsets,
prevalences,
directionUpper,
piTreatments,
piControls,
stratifiedAnalysis,
plannedSubjects,
typeOfSelection,
effectMeasure,
adaptations,
epsilonValue,
rValue,
threshold,
allocationRatioPlanned,
minNumberOfSubjectsPerStage,
maxNumberOfSubjectsPerStage,
conditionalPower,
piTreatmentH1,
piControlH1,
calcSubjectsFunction,
calcSubjectsFunctionIsUserDefined,
selectPopulationsFunction) {
kMax <- length(plannedSubjects)
pMax <- length(piTreatments)
gMax <- log(length(piTreatments), 2) + 1
subjectsPerStage <- matrix(NA_real_, nrow = pMax, ncol = kMax)
simEventsTreatment <- matrix(NA_real_, nrow = pMax, ncol = kMax)
simEventsControl <- matrix(NA_real_, nrow = pMax, ncol = kMax)
populationSubjectsPerStage <- matrix(NA_real_, nrow = gMax, ncol = kMax)
testStatistics <- matrix(NA_real_, nrow = gMax, ncol = kMax)
overallTestStatistics <- matrix(NA_real_, nrow = gMax, ncol = kMax)
separatePValues <- matrix(NA_real_, nrow = gMax, ncol = kMax)
conditionalCriticalValue <- rep(NA_real_, kMax - 1)
conditionalPowerPerStage <- rep(NA_real_, kMax)
selectedPopulations <- matrix(FALSE, nrow = gMax, ncol = kMax)
selectedSubsets <- matrix(FALSE, nrow = pMax, ncol = kMax)
selectedPopulations[, 1] <- TRUE
selectedSubsets[, 1] <- TRUE
adjustedPValues <- rep(NA_real_, kMax)
overallRatesTreatment <- matrix(NA_real_, nrow = gMax, ncol = kMax)
overallRatesControl <- matrix(NA_real_, nrow = gMax, ncol = kMax)
overallEffectSizes <- matrix(NA_real_, nrow = gMax, ncol = kMax)
if (.isTrialDesignFisher(design)) {
weights <- .getWeightsFisher(design)
} else if (.isTrialDesignInverseNormal(design)) {
weights <- .getWeightsInverseNormal(design)
}
for (k in 1:kMax) {
const <- allocationRatioPlanned[k]
selectedSubsets[, k] <- .createSelectedSubsets(k, selectedPopulations)
if (k == 1) {
# subjectsPerStage[, k] <- stats::rmultinom(1, plannedSubjects[k], prevalences)
subjectsPerStage[, k] <- plannedSubjects[k] * prevalences
} else {
prevSelected <- prevalences / sum(prevalences[selectedSubsets[, k]])
prevSelected[!selectedSubsets[, k]] <- 0
if (sum(prevSelected, na.rm = TRUE) > 0) {
# subjectsPerStage[, k] <- stats::rmultinom(1, plannedSubjects[k] - plannedSubjects[k - 1], prevSelected)
subjectsPerStage[, k] <- (plannedSubjects[k] - plannedSubjects[k - 1]) * prevSelected
} else {
break
}
}
selsubs <- !is.na(subjectsPerStage[, k]) & subjectsPerStage[, k] > 0
if (any(round(subjectsPerStage[selsubs, k] * const / (1 + const)) < 1) ||
any(round(subjectsPerStage[selsubs, k] / (1 + const)) < 1)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"at least one sample size specification too small to create simulation results, ",
"e.g., due to small prevalences of subsets"
)
}
simEventsTreatment[selsubs, k] <- stats::rbinom(
rep(1, sum(selsubs)),
round(subjectsPerStage[selsubs, k] * const / (1 + const)), piTreatments[selsubs]
)
simEventsControl[selsubs, k] <- stats::rbinom(
rep(1, sum(selsubs)),
round(subjectsPerStage[selsubs, k] / (1 + const)), piControls[selsubs]
)
if (gMax == 1) {
rm <- (simEventsControl[1, k] + simEventsTreatment[1, k]) / subjectsPerStage[1, k]
if (rm <= 0 || rm >= 1) {
testStatistics[1, k] <- 0
} else {
testStatistics[1, k] <- (2 * directionUpper - 1) *
(simEventsTreatment[1, k] * (1 + const) / const - simEventsControl[1, k] *
(1 + const)) / subjectsPerStage[1, k] /
sqrt(rm * (1 - rm)) * sqrt(subjectsPerStage[1, k] * const / (1 + const)^2)
}
populationSubjectsPerStage[1, k] <- subjectsPerStage[1, k]
overallRatesTreatment[1, k] <- sum(simEventsTreatment[1, 1:k]) /
round(const / (1 + const) * sum(subjectsPerStage[1, 1:k]))
overallRatesControl[1, k] <- sum(simEventsControl[1, 1:k]) /
round(1 / (1 + const) * sum(subjectsPerStage[1, 1:k]))
overallEffectSizes[1, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[1, k] - overallRatesControl[1, k])
rm <- sum(simEventsControl[1, 1:k] + simEventsTreatment[1, 1:k]) / sum(subjectsPerStage[1, 1:k])
if (rm == 0 || rm == 1) {
overallTestStatistics[1, k] <- 0
} else {
overallTestStatistics[1, k] <- overallEffectSizes[1, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[1, 1:k]) * const / (1 + const)^2)
}
} else if (gMax == 2) {
# Population S1
rm <- (simEventsControl[1, k] + simEventsTreatment[1, k]) / subjectsPerStage[1, k]
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[1, k] <- 0
} else {
testStatistics[1, k] <- (2 * directionUpper - 1) *
(simEventsTreatment[1, k] * (1 + const) / const - simEventsControl[1, k] *
(1 + const)) / subjectsPerStage[1, k] /
sqrt(rm * (1 - rm)) * sqrt(subjectsPerStage[1, k] * const / (1 + const)^2)
}
}
populationSubjectsPerStage[1, k] <- subjectsPerStage[1, k]
overallRatesTreatment[1, k] <- sum(simEventsTreatment[1, 1:k]) /
round(const / (1 + const) * sum(subjectsPerStage[1, 1:k]))
overallRatesControl[1, k] <- sum(simEventsControl[1, 1:k]) /
round(1 / (1 + const) * sum(subjectsPerStage[1, 1:k]))
overallEffectSizes[1, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[1, k] - overallRatesControl[1, k])
rm <- sum(simEventsControl[1, 1:k] + simEventsTreatment[1, 1:k]) / sum(subjectsPerStage[1, 1:k])
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[1, k] <- 0
} else {
overallTestStatistics[1, k] <- overallEffectSizes[1, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[1, 1:k]) * const / (1 + const)^2)
}
}
# Full population
if (stratifiedAnalysis) {
rm <- (simEventsControl[1:2, k] + simEventsTreatment[1:2, k]) / subjectsPerStage[1:2, k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[2, k] <- 0
} else {
testStatistics[2, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[1:2, k] * (simEventsTreatment[1:2, k] *
(1 + const) / const - simEventsControl[1:2, k] * (1 + const)) /
subjectsPerStage[1:2, k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[1:2, k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[1:2, k] + simEventsTreatment[1:2, k], na.rm = TRUE) /
sum(subjectsPerStage[1:2, k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[2, k] <- 0
} else {
testStatistics[2, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[1:2, k] * (1 + const) / const - simEventsControl[1:2, k] *
(1 + const), na.rm = TRUE) /
sum(subjectsPerStage[1:2, k], na.rm = TRUE) /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[1:2, k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[2, k] <- sum(subjectsPerStage[1:2, k], na.rm = TRUE)
overallRatesTreatment[2, k] <- sum(simEventsTreatment[1:2, 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[1:2, 1:k], na.rm = TRUE))
overallRatesControl[2, k] <- sum(simEventsControl[1:2, 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[1:2, 1:k], na.rm = TRUE))
overallEffectSizes[2, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[2, k] - overallRatesControl[2, k])
rm <- sum(simEventsControl[1:2, 1:k] + simEventsTreatment[1:2, 1:k], na.rm = TRUE) /
sum(subjectsPerStage[1:2, 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[2, k] <- 0
} else {
overallTestStatistics[2, k] <- overallEffectSizes[2, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[1:2, 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
} else if (gMax == 3) {
# Population S1
if (stratifiedAnalysis) {
rm <- (simEventsControl[c(1, 3), k] + simEventsTreatment[c(1, 3), k]) / subjectsPerStage[c(1, 3), k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[1, k] <- 0
} else {
testStatistics[1, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[c(1, 3), k] * (simEventsTreatment[c(1, 3), k] *
(1 + const) / const - simEventsControl[c(1, 3), k] * (1 + const)) /
subjectsPerStage[c(1, 3), k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[c(1, 3), k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[c(1, 3), k] + simEventsTreatment[c(1, 3), k], na.rm = TRUE) /
sum(subjectsPerStage[c(1, 3), k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[1, k] <- 0
} else {
testStatistics[1, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[c(1, 3), k] * (1 + const) / const - simEventsControl[c(1, 3), k] *
(1 + const), na.rm = TRUE) /
sum(subjectsPerStage[c(1, 3), k], na.rm = TRUE) /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(1, 3), k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[1, k] <- sum(subjectsPerStage[c(1, 3), k], na.rm = TRUE)
overallRatesTreatment[1, k] <- sum(simEventsTreatment[c(1, 3), 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[c(1, 3), 1:k], na.rm = TRUE))
overallRatesControl[1, k] <- sum(simEventsControl[c(1, 3), 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[c(1, 3), 1:k], na.rm = TRUE))
overallEffectSizes[1, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[1, k] - overallRatesControl[1, k])
rm <- sum(simEventsControl[c(1, 3), 1:k] + simEventsTreatment[c(1, 3), 1:k], na.rm = TRUE) /
sum(subjectsPerStage[c(1, 3), 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[1, k] <- 0
} else {
overallTestStatistics[1, k] <- overallEffectSizes[1, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(1, 3), 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
# Population S2
if (stratifiedAnalysis) {
rm <- (simEventsControl[c(2, 3), k] + simEventsTreatment[c(2, 3), k]) / subjectsPerStage[c(2, 3), k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[2, k] <- 0
} else {
testStatistics[2, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[c(2, 3), k] * (simEventsTreatment[c(2, 3), k] *
(1 + const) / const - simEventsControl[c(2, 3), k] * (1 + const)) /
subjectsPerStage[c(2, 3), k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[c(2, 3), k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[c(2, 3), k] + simEventsTreatment[c(2, 3), k], na.rm = TRUE) /
sum(subjectsPerStage[c(2, 3), k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[2, k] <- 0
} else {
testStatistics[2, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[c(2, 3), k] * (1 + const) / const - simEventsControl[c(2, 3), k] *
(1 + const), na.rm = TRUE) /
sum(subjectsPerStage[c(2, 3), k], na.rm = TRUE) /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(2, 3), k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[2, k] <- sum(subjectsPerStage[c(2, 3), k], na.rm = TRUE)
overallRatesTreatment[2, k] <- sum(simEventsTreatment[c(2, 3), 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[c(2, 3), 1:k], na.rm = TRUE))
overallRatesControl[2, k] <- sum(simEventsControl[c(2, 3), 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[c(2, 3), 1:k], na.rm = TRUE))
overallEffectSizes[2, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[2, k] - overallRatesControl[2, k])
rm <- sum(simEventsControl[c(2, 3), 1:k] + simEventsTreatment[c(2, 3), 1:k], na.rm = TRUE) /
sum(subjectsPerStage[c(2, 3), 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[2, k] <- 0
} else {
overallTestStatistics[2, k] <- overallEffectSizes[2, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(2, 3), 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
# Full population
if (stratifiedAnalysis) {
rm <- (simEventsControl[1:4, k] + simEventsTreatment[1:4, k]) / subjectsPerStage[1:4, k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[3, k] <- 0
} else {
testStatistics[3, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[1:4, k] * (simEventsTreatment[1:4, k] *
(1 + const) / const - simEventsControl[1:4, k] * (1 + const)) /
subjectsPerStage[1:4, k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[1:4, k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[1:4, k] + simEventsTreatment[1:4, k], na.rm = TRUE) /
sum(subjectsPerStage[1:4, k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[3, k] <- 0
} else {
testStatistics[3, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[1:4, k] * (1 + const) / const - simEventsControl[1:4, k] *
(1 + const), na.rm = TRUE) /
sum(subjectsPerStage[1:4, k], na.rm = TRUE) /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[1:4, k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[3, k] <- sum(subjectsPerStage[1:4, k], na.rm = TRUE)
overallRatesTreatment[3, k] <- sum(simEventsTreatment[1:4, 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[1:4, 1:k], na.rm = TRUE))
overallRatesControl[3, k] <- sum(simEventsControl[1:4, 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[1:4, 1:k], na.rm = TRUE))
overallEffectSizes[3, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[3, k] - overallRatesControl[3, k])
rm <- sum(simEventsControl[1:4, 1:k] + simEventsTreatment[1:4, 1:k], na.rm = TRUE) /
sum(subjectsPerStage[1:4, 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[3, k] <- 0
} else {
overallTestStatistics[3, k] <- overallEffectSizes[3, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[1:4, 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
} else if (gMax == 4) {
# Population S1
if (stratifiedAnalysis) {
rm <- (simEventsControl[c(1, 4, 5, 7), k] + simEventsTreatment[c(1, 4, 5, 7), k]) / subjectsPerStage[c(1, 4, 5, 7), k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[1, k] <- 0
} else {
testStatistics[1, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[c(1, 4, 5, 7), k] * (simEventsTreatment[c(1, 4, 5, 7), k] *
(1 + const) / const - simEventsControl[c(1, 4, 5, 7), k] * (1 + const)) /
subjectsPerStage[c(1, 4, 5, 7), k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[c(1, 4, 5, 7), k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[c(1, 4, 5, 7), k] + simEventsTreatment[c(1, 4, 5, 7), k], na.rm = TRUE) /
sum(subjectsPerStage[c(1, 4, 5, 7), k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[1, k] <- 0
} else {
testStatistics[1, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[c(1, 4, 5, 7), k] * (1 + const) / const -
simEventsControl[c(1, 4, 5, 7), k] * (1 + const), na.rm = TRUE) /
sum(subjectsPerStage[c(1, 4, 5, 7), k], na.rm = TRUE) /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(1, 4, 5, 7), k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[1, k] <- sum(subjectsPerStage[c(1, 4, 5, 7), k], na.rm = TRUE)
overallRatesTreatment[1, k] <- sum(simEventsTreatment[c(1, 4, 5, 7), 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[c(1, 4, 5, 7), 1:k], na.rm = TRUE))
overallRatesControl[1, k] <- sum(simEventsControl[c(1, 4, 5, 7), 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[c(1, 4, 5, 7), 1:k], na.rm = TRUE))
overallEffectSizes[1, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[1, k] - overallRatesControl[1, k])
rm <- sum(simEventsControl[c(1, 4, 5, 7), 1:k] + simEventsTreatment[c(1, 4, 5, 7), 1:k], na.rm = TRUE) /
sum(subjectsPerStage[c(1, 4, 5, 7), 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[1, k] <- 0
} else {
overallTestStatistics[1, k] <- overallEffectSizes[1, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(1, 4, 5, 7), 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
# Population S2
if (stratifiedAnalysis) {
rm <- (simEventsControl[c(2, 4, 6, 7), k] + simEventsTreatment[c(2, 4, 6, 7), k]) / subjectsPerStage[c(2, 4, 6, 7), k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[2, k] <- 0
} else {
testStatistics[2, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[c(2, 4, 6, 7), k] * (simEventsTreatment[c(2, 4, 6, 7), k] *
(1 + const) / const - simEventsControl[c(2, 4, 6, 7), k] * (1 + const)) /
subjectsPerStage[c(2, 4, 6, 7), k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[c(2, 4, 6, 7), k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[c(2, 4, 6, 7), k] + simEventsTreatment[c(2, 4, 6, 7), k], na.rm = TRUE) /
sum(subjectsPerStage[c(2, 4, 6, 7), k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[2, k] <- 0
} else {
testStatistics[2, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[c(2, 4, 6, 7), k] * (1 + const) / const - simEventsControl[c(2, 4, 6, 7), k] *
(1 + const), na.rm = TRUE) /
sum(subjectsPerStage[c(2, 4, 6, 7), k], na.rm = TRUE) / sqrt(rm * (1 - rm)) *
sqrt(sum(subjectsPerStage[c(2, 4, 6, 7), k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[2, k] <- sum(subjectsPerStage[c(2, 4, 6, 7), k], na.rm = TRUE)
overallRatesTreatment[2, k] <- sum(simEventsTreatment[c(2, 4, 6, 7), 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[c(2, 4, 6, 7), 1:k], na.rm = TRUE))
overallRatesControl[2, k] <- sum(simEventsControl[c(2, 4, 6, 7), 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[c(2, 4, 6, 7), 1:k], na.rm = TRUE))
overallEffectSizes[2, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[2, k] - overallRatesControl[2, k])
rm <- sum(simEventsControl[c(2, 4, 6, 7), 1:k] + simEventsTreatment[c(2, 4, 6, 7), 1:k], na.rm = TRUE) /
sum(subjectsPerStage[c(2, 4, 6, 7), 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[2, k] <- 0
} else {
overallTestStatistics[2, k] <- overallEffectSizes[2, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(2, 4, 6, 7), 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
# Population S3
if (stratifiedAnalysis) {
rm <- (simEventsControl[c(3, 5, 6, 7), k] + simEventsTreatment[c(3, 5, 6, 7), k]) / subjectsPerStage[c(3, 5, 6, 7), k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[3, k] <- 0
} else {
testStatistics[3, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[c(3, 5, 6, 7), k] * (simEventsTreatment[c(3, 5, 6, 7), k] *
(1 + const) / const - simEventsControl[c(3, 5, 6, 7), k] * (1 + const)) /
subjectsPerStage[c(3, 5, 6, 7), k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[c(3, 5, 6, 7), k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[c(3, 5, 6, 7), k] + simEventsTreatment[c(3, 5, 6, 7), k], na.rm = TRUE) /
sum(subjectsPerStage[c(3, 5, 6, 7), k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[3, k] <- 0
} else {
testStatistics[3, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[c(3, 5, 6, 7), k] * (1 + const) / const -
simEventsControl[c(3, 5, 6, 7), k] * (1 + const), na.rm = TRUE) /
sum(subjectsPerStage[c(3, 5, 6, 7), k], na.rm = TRUE) / sqrt(rm * (1 - rm)) *
sqrt(sum(subjectsPerStage[c(3, 5, 6, 7), k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[3, k] <- sum(subjectsPerStage[c(3, 5, 6, 7), k], na.rm = TRUE)
overallRatesTreatment[3, k] <- sum(simEventsTreatment[c(3, 5, 6, 7), 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[c(3, 5, 6, 7), 1:k], na.rm = TRUE))
overallRatesControl[3, k] <- sum(simEventsControl[c(3, 5, 6, 7), 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[c(3, 5, 6, 7), 1:k], na.rm = TRUE))
overallEffectSizes[3, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[3, k] - overallRatesControl[3, k])
rm <- sum(simEventsControl[c(3, 5, 6, 7), 1:k] + simEventsTreatment[c(3, 5, 6, 7), 1:k], na.rm = TRUE) /
sum(subjectsPerStage[c(3, 5, 6, 7), 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[3, k] <- 0
} else {
overallTestStatistics[3, k] <- overallEffectSizes[3, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[c(3, 5, 6, 7), 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
# Full population
if (stratifiedAnalysis) {
rm <- (simEventsControl[1:8, k] + simEventsTreatment[1:8, k]) / subjectsPerStage[1:8, k]
rm[!is.na(rm) & (rm <= 0 | rm >= 1)] <- 0
if (!all(is.na(rm))) {
if (all(na.omit(rm) == 0) || all(na.omit(rm) == 1)) {
testStatistics[4, k] <- 0
} else {
testStatistics[4, k] <- sqrt(const) / (1 + const) * (2 * directionUpper - 1) *
sum(subjectsPerStage[1:8, k] * (simEventsTreatment[1:8, k] *
(1 + const) / const - simEventsControl[1:8, k] * (1 + const)) /
subjectsPerStage[1:8, k], na.rm = TRUE) / sqrt(sum(rm * (1 - rm) *
subjectsPerStage[1:8, k], na.rm = TRUE))
}
}
} else {
rm <- sum(simEventsControl[1:8, k] + simEventsTreatment[1:8, k], na.rm = TRUE) /
sum(subjectsPerStage[1:8, k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
testStatistics[4, k] <- 0
} else {
testStatistics[4, k] <- (2 * directionUpper - 1) *
sum(simEventsTreatment[1:8, k] * (1 + const) / const - simEventsControl[1:8, k] *
(1 + const), na.rm = TRUE) /
sum(subjectsPerStage[1:8, k], na.rm = TRUE) / sqrt(rm * (1 - rm)) *
sqrt(sum(subjectsPerStage[1:8, k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
populationSubjectsPerStage[4, k] <- sum(subjectsPerStage[1:8, k], na.rm = TRUE)
overallRatesTreatment[4, k] <- sum(simEventsTreatment[1:8, 1:k], na.rm = TRUE) /
round(const / (1 + const) * sum(subjectsPerStage[1:8, 1:k], na.rm = TRUE))
overallRatesControl[4, k] <- sum(simEventsControl[1:8, 1:k], na.rm = TRUE) /
round(1 / (1 + const) * sum(subjectsPerStage[1:8, 1:k], na.rm = TRUE))
overallEffectSizes[4, k] <-
(2 * directionUpper - 1) * (overallRatesTreatment[4, k] - overallRatesControl[4, k])
rm <- sum(simEventsControl[1:8, 1:k] + simEventsTreatment[1:8, 1:k], na.rm = TRUE) /
sum(subjectsPerStage[1:8, 1:k], na.rm = TRUE)
if (!is.na(rm)) {
if (rm <= 0 || rm >= 1) {
overallTestStatistics[4, k] <- 0
} else {
overallTestStatistics[4, k] <- overallEffectSizes[4, k] /
sqrt(rm * (1 - rm)) * sqrt(sum(subjectsPerStage[1:8, 1:k], na.rm = TRUE) * const / (1 + const)^2)
}
}
}
testStatistics[!selectedPopulations[, k], k] <- NA_real_
overallEffectSizes[!selectedPopulations[, k], k] <- NA_real_
overallTestStatistics[!selectedPopulations[, k], k] <- NA_real_
separatePValues[, k] <- 1 - stats::pnorm(testStatistics[, k])
if (k < kMax) {
if (colSums(selectedPopulations)[k] == 0) {
break
}
# Bonferroni adjustment
adjustedPValues[k] <- min(min(separatePValues[, k], na.rm = TRUE) *
colSums(selectedPopulations)[k], 1 - 1e-12)
# conditional critical value to reject the null hypotheses at the next stage of the trial
if (.isTrialDesignFisher(design)) {
conditionalCriticalValue[k] <- .getOneMinusQNorm(min((design$criticalValues[k + 1] /
prod(adjustedPValues[1:k]^weights[1:k]))^(1 / weights[k + 1]), 1 - 1e-07))
} else {
if (design$criticalValues[k + 1] >= 6) {
conditionalCriticalValue[k] <- Inf
} else {
conditionalCriticalValue[k] <- (design$criticalValues[k + 1] * sqrt(design$informationRates[k + 1]) -
.getOneMinusQNorm(adjustedPValues[1:k]) %*% weights[1:k]) /
sqrt(design$informationRates[k + 1] - design$informationRates[k])
}
}
if (adaptations[k]) {
if (effectMeasure == "testStatistic") {
selectedPopulations[, k + 1] <- (selectedPopulations[, k] &
.selectPopulations(
k, overallTestStatistics[, k] + runif(gMax, -1e-05, 1e-05),
typeOfSelection, epsilonValue, rValue, threshold, selectPopulationsFunction
))
} else if (effectMeasure == "effectEstimate") {
selectedPopulations[, k + 1] <- (selectedPopulations[, k] &
.selectPopulations(
k, overallEffectSizes[, k] + runif(gMax, -1e-05, 1e-05),
typeOfSelection, epsilonValue, rValue, threshold, selectPopulationsFunction
))
}
newSubjects <- calcSubjectsFunction(
stage = k + 1, # to be consistent with non-enrichment situation, cf. line 40
directionUpper = directionUpper,
conditionalPower = conditionalPower,
conditionalCriticalValue = conditionalCriticalValue,
plannedSubjects = plannedSubjects,
allocationRatioPlanned = allocationRatioPlanned,
selectedPopulations = selectedPopulations,
piTreatmentH1 = piTreatmentH1,
piControlH1 = piControlH1,
overallRatesTreatment = overallRatesTreatment,
overallRatesControl = overallRatesControl,
minNumberOfSubjectsPerStage = minNumberOfSubjectsPerStage,
maxNumberOfSubjectsPerStage = maxNumberOfSubjectsPerStage
)
if (is.null(newSubjects) || length(newSubjects) != 1 || !is.numeric(newSubjects) || is.na(newSubjects)) {
stop(
C_EXCEPTION_TYPE_ILLEGAL_ARGUMENT,
"'calcSubjectsFunction' returned an illegal or undefined result (", newSubjects, "); ",
"the output must be a single numeric value"
)
}
if (!is.na(conditionalPower) || calcSubjectsFunctionIsUserDefined) {
plannedSubjects[(k + 1):kMax] <- plannedSubjects[k] + cumsum(rep(newSubjects, kMax - k))
}
} else {
selectedPopulations[, k + 1] <- selectedPopulations[, k]
}
if (is.na(piControlH1)) {
pi2H1 <- overallRatesControl[, k]
} else {
pi2H1 <- piControlH1
}
if (is.na(piTreatmentH1)) {
pi1H1 <- overallRatesTreatment[, k]
} else {
pi1H1 <- piTreatmentH1
}
pim <- (allocationRatioPlanned[k] * pi1H1 + pi2H1) / (1 + allocationRatioPlanned[k])
if (any(pi1H1 * (1 - pi1H1) + pi2H1 * (1 - pi2H1) == 0)) {
thetaStandardized <- 0
} else {
thetaStandardized <- sqrt(allocationRatioPlanned[k]) / (1 + allocationRatioPlanned[k]) * (
(pi1H1 - pi2H1) * sqrt(1 + allocationRatioPlanned[k]) /
sqrt(pi1H1 * (1 - pi1H1) + allocationRatioPlanned[k] * pi2H1 * (1 - pi2H1)) +
sign(pi1H1 - pi2H1) * conditionalCriticalValue[k] *
(1 - sqrt(pim * (1 - pim) + allocationRatioPlanned[k] * pim * (1 - pim)) /
sqrt(pi1H1 * (1 - pi1H1) + allocationRatioPlanned[k] * pi2H1 * (1 - pi2H1))) *
(1 + allocationRatioPlanned[k]) / sqrt(allocationRatioPlanned[k] * (plannedSubjects[k + 1] - plannedSubjects[k]))
)
}
thetaStandardized <- (2 * directionUpper - 1) * thetaStandardized
thetaStandardized <- min(thetaStandardized, na.rm = TRUE)
conditionalPowerPerStage[k] <- 1 - stats::pnorm(conditionalCriticalValue[k] -
thetaStandardized * sqrt(plannedSubjects[k + 1] - plannedSubjects[k]))
}
}
return(list(
subjectsPerStage = subjectsPerStage,
populationSubjectsPerStage = populationSubjectsPerStage,
allocationRatioPlanned = allocationRatioPlanned,
overallEffectSizes = overallEffectSizes,
testStatistics = testStatistics,
directionUpper = directionUpper,
overallTestStatistics = overallTestStatistics,
overallRatesControl = overallRatesControl,
overallRatesTreatment = overallRatesTreatment,
separatePValues = separatePValues,
conditionalCriticalValue = conditionalCriticalValue,
conditionalPowerPerStage = conditionalPowerPerStage,
selectedPopulations = selectedPopulations
))
}
#'
#' @title
#' Get Simulation Enrichment Rates
#'
#' @description
#' Returns the simulated power, stopping and selection probabilities, conditional power,
#' and expected sample size for testing rates in an enrichment design testing situation.
#'
#' @param piControlH1 If specified, the assumed probabilities in the control arm
#' under which the sample size recalculation was performed
#' and the conditional power was calculated.
#' @param piTreatmentH1 If specified, the assumed probabilities in the active arm
#' under which the sample size recalculation was performed
#' and the conditional power was calculated.
#' @inheritParams param_intersectionTest_Enrichment
#' @inheritParams param_typeOfSelection
#' @inheritParams param_effectMeasure
#' @inheritParams param_adaptations
#' @inheritParams param_threshold
#' @inheritParams param_effectList
#' @inheritParams param_successCriterion
#' @inheritParams param_typeOfSelection
#' @inheritParams param_design_with_default
#' @inheritParams param_directionUpper
#' @inheritParams param_allocationRatioPlanned
#' @inheritParams param_plannedSubjects
#' @inheritParams param_minNumberOfSubjectsPerStage
#' @inheritParams param_maxNumberOfSubjectsPerStage
#' @inheritParams param_conditionalPowerSimulation
#' @inheritParams param_maxNumberOfIterations
#' @inheritParams param_calcSubjectsFunction
#' @inheritParams param_selectPopulationsFunction
#' @inheritParams param_rValue
#' @inheritParams param_epsilonValue
#' @inheritParams param_seed
#' @inheritParams param_three_dots
#' @inheritParams param_showStatistics
#' @inheritParams param_stratifiedAnalysis
#'
#' @details
#' At given design the function simulates the power, stopping probabilities,
#' selection probabilities, and expected sample size at given number of subjects,
#' parameter configuration, and treatment arm selection rule in the enrichment situation.
#' An allocation ratio can be specified referring to the ratio of number of
#' subjects in the active treatment groups as compared to the control group.
#'
#' The definition of \code{piTreatmentH1} and/or \code{piControlH1} makes only sense if \code{kMax} > 1
#' and if \code{conditionalPower}, \code{minNumberOfSubjectsPerStage}, and
#' \code{maxNumberOfSubjectsPerStage} (or \code{calcSubjectsFunction}) are defined.
#'
#' \code{calcSubjectsFunction}\cr
#' This function returns the number of subjects at given conditional power and
#' conditional critical value for specified testing situation.
#' The function might depend on the variables
#' \code{stage},
#' \code{selectedPopulations},
#' \code{directionUpper},
#' \code{plannedSubjects},
#' \code{allocationRatioPlanned},
#' \code{minNumberOfSubjectsPerStage},
#' \code{maxNumberOfSubjectsPerStage},
#' \code{conditionalPower},
#' \code{conditionalCriticalValue},
#' \code{overallRatesTreatment},
#' \code{overallRatesControl},
#' \code{piTreatmentH1}, and
#' \code{piControlH1}.
#' The function has to contain the three-dots argument '...' (see examples).
#'
#' @template return_object_simulation_results
#' @template how_to_get_help_for_generics
#'
#' @template examples_get_simulation_enrichment_rates
#'
#' @export
#'
getSimulationEnrichmentRates <- function(design = NULL, ...,
effectList = NULL,
intersectionTest = c("Simes", "SpiessensDebois", "Bonferroni", "Sidak"), # C_INTERSECTION_TEST_ENRICHMENT_DEFAULT
stratifiedAnalysis = TRUE, # C_STRATIFIED_ANALYSIS_DEFAULT,
directionUpper = TRUE, # C_DIRECTION_UPPER_DEFAULT
adaptations = NA,
typeOfSelection = c("best", "rBest", "epsilon", "all", "userDefined"), # C_TYPE_OF_SELECTION_DEFAULT
effectMeasure = c("effectEstimate", "testStatistic"), # C_EFFECT_MEASURE_DEFAULT
successCriterion = c("all", "atLeastOne"), # C_SUCCESS_CRITERION_DEFAULT
epsilonValue = NA_real_,
rValue = NA_real_,
threshold = -Inf,
plannedSubjects = NA_real_,
allocationRatioPlanned = NA_real_,
minNumberOfSubjectsPerStage = NA_real_,
maxNumberOfSubjectsPerStage = NA_real_,
conditionalPower = NA_real_,
piTreatmentH1 = NA_real_,
piControlH1 = NA_real_,
maxNumberOfIterations = 1000L, # C_MAX_SIMULATION_ITERATIONS_DEFAULT
seed = NA_real_,
calcSubjectsFunction = NULL,
selectPopulationsFunction = NULL,
showStatistics = FALSE) {
if (is.null(design)) {
design <- .getDefaultDesign(..., type = "simulation")
.warnInCaseOfUnknownArguments(
functionName = "getSimulationEnrichmentRates",
ignore = c(.getDesignArgumentsToIgnoreAtUnknownArgumentCheck(
design,
powerCalculationEnabled = TRUE
), "showStatistics"), ...
)
} else {
.assertIsTrialDesignInverseNormalOrFisher(design)
.warnInCaseOfUnknownArguments(functionName = "getSimulationEnrichmentRates", ignore = "showStatistics", ...)
.warnInCaseOfTwoSidedPowerArgument(...)
}
.assertIsOneSidedDesign(design, designType = "enrichment", engineType = "simulation")
calcSubjectsFunctionIsUserDefined <- !is.null(calcSubjectsFunction)
simulationResults <- .createSimulationResultsEnrichmentObject(
design = design,
effectList = effectList,
intersectionTest = intersectionTest,
stratifiedAnalysis = stratifiedAnalysis,
directionUpper = directionUpper, # rates + survival only
adaptations = adaptations,
typeOfSelection = typeOfSelection,
effectMeasure = effectMeasure,
successCriterion = successCriterion,
epsilonValue = epsilonValue,
rValue = rValue,
threshold = threshold,
plannedSubjects = plannedSubjects, # means + rates only
allocationRatioPlanned = allocationRatioPlanned,
minNumberOfSubjectsPerStage = minNumberOfSubjectsPerStage, # means + rates only
maxNumberOfSubjectsPerStage = maxNumberOfSubjectsPerStage, # means + rates only
conditionalPower = conditionalPower,
piTreatmentH1 = piTreatmentH1, # rates only
piControlH1 = piControlH1, # rates only
maxNumberOfIterations = maxNumberOfIterations,
seed = seed,
calcSubjectsFunction = calcSubjectsFunction, # means + rates only
selectPopulationsFunction = selectPopulationsFunction,
showStatistics = showStatistics,
endpoint = "rates"
)
design <- simulationResults$.design
successCriterion <- simulationResults$successCriterion
effectMeasure <- simulationResults$effectMeasure
adaptations <- simulationResults$adaptations
gMax <- simulationResults$populations
kMax <- simulationResults$.design$kMax
intersectionTest <- simulationResults$intersectionTest
typeOfSelection <- simulationResults$typeOfSelection
effectList <- simulationResults$effectList
piTreatmentH1 <- simulationResults$piTreatmentH1 # rates only
piControlH1 <- simulationResults$piControlH1 # rates only
conditionalPower <- simulationResults$conditionalPower
minNumberOfSubjectsPerStage <- simulationResults$minNumberOfSubjectsPerStage
maxNumberOfSubjectsPerStage <- simulationResults$maxNumberOfSubjectsPerStage
allocationRatioPlanned <- simulationResults$allocationRatioPlanned
calcSubjectsFunction <- simulationResults$calcSubjectsFunction
if (length(allocationRatioPlanned) == 1) {
allocationRatioPlanned <- rep(allocationRatioPlanned, kMax)
}
indices <- .getIndicesOfClosedHypothesesSystemForSimulation(gMax = gMax)
cols <- nrow(effectList$piTreatments)
simulatedSelections <- array(0, dim = c(kMax, cols, gMax))
simulatedRejections <- array(0, dim = c(kMax, cols, gMax))
simulatedNumberOfPopulations <- matrix(0, nrow = kMax, ncol = cols)
simulatedSubjectsPerStage <- array(0, dim = c(kMax, cols, 2^(gMax - 1)))
simulatedSuccessStopping <- matrix(0, nrow = kMax, ncol = cols)
simulatedFutilityStopping <- matrix(0, nrow = kMax - 1, ncol = cols)
simulatedConditionalPower <- matrix(0, nrow = kMax, ncol = cols)
simulatedRejectAtLeastOne <- rep(0, cols)
expectedNumberOfSubjects <- rep(0, cols)
iterations <- matrix(0, nrow = kMax, ncol = cols)
len <- maxNumberOfIterations * kMax * gMax * cols
dataIterationNumber <- rep(NA_real_, len)
dataStageNumber <- rep(NA_real_, len)
dataPopulationNumber <- rep(NA_real_, len)
dataEffect <- rep(NA_real_, len)
dataSubjectsPopulation <- rep(NA_real_, len)
dataSubjectsActivePopulation <- rep(NA_real_, len)
dataNumberOfSubjects <- rep(NA_real_, len)
dataNumberOfCumulatedSubjects <- rep(NA_real_, len)
dataRejectPerStage <- rep(NA, len)
dataFutilityStop <- rep(NA_real_, len)
dataSuccessStop <- rep(NA, len)
dataFutilityStop <- rep(NA, len)
dataTestStatistics <- rep(NA_real_, len)
dataConditionalCriticalValue <- rep(NA_real_, len)
dataConditionalPowerAchieved <- rep(NA_real_, len)
dataEffectEstimate <- rep(NA_real_, len)
dataPValuesSeparate <- rep(NA_real_, len)
piControls <- effectList$piControls
if (length(piControls) == 1) {
piControls <- rep(piControls, ncol(effectList$piTreatments))
}
index <- 1
for (i in 1:cols) {
for (j in 1:maxNumberOfIterations) {
stageResults <- .getSimulatedStageRatesEnrichment(
design = design,
subsets = effectList$subsets,
prevalences = effectList$prevalences,
piTreatments = effectList$piTreatments[i, ],
piControls = piControls,
directionUpper = directionUpper,
stratifiedAnalysis = stratifiedAnalysis,
plannedSubjects = plannedSubjects,
typeOfSelection = typeOfSelection,
effectMeasure = effectMeasure,
adaptations = adaptations,
epsilonValue = epsilonValue,
rValue = rValue,
threshold = threshold,
allocationRatioPlanned = allocationRatioPlanned,
minNumberOfSubjectsPerStage = minNumberOfSubjectsPerStage,
maxNumberOfSubjectsPerStage = maxNumberOfSubjectsPerStage,
conditionalPower = conditionalPower,
piTreatmentH1 = piTreatmentH1,
piControlH1 = piControlH1,
calcSubjectsFunction = calcSubjectsFunction,
calcSubjectsFunctionIsUserDefined = calcSubjectsFunctionIsUserDefined,
selectPopulationsFunction = selectPopulationsFunction
)
closedTest <- .performClosedCombinationTestForSimulationEnrichment(
stageResults = stageResults,
design = design, indices = indices,
intersectionTest = intersectionTest, successCriterion = successCriterion
)
rejectAtSomeStage <- FALSE
rejectedPopulationsBefore <- rep(FALSE, gMax)
for (k in 1:kMax) {
simulatedRejections[k, i, ] <- simulatedRejections[k, i, ] +
(closedTest$rejected[, k] & closedTest$selectedPopulations[1:gMax, k] | rejectedPopulationsBefore)
simulatedSelections[k, i, ] <- simulatedSelections[k, i, ] + closedTest$selectedPopulations[, k]
simulatedNumberOfPopulations[k, i] <- simulatedNumberOfPopulations[k, i] +
sum(closedTest$selectedPopulations[, k])
if (!any(is.na(closedTest$successStop))) {
simulatedSuccessStopping[k, i] <- simulatedSuccessStopping[k, i] + closedTest$successStop[k]
}
if ((kMax > 1) && (k < kMax)) {
if (!any(is.na(closedTest$futilityStop))) {
simulatedFutilityStopping[k, i] <- simulatedFutilityStopping[k, i] +
(closedTest$futilityStop[k] && !closedTest$successStop[k])
}
if (!closedTest$successStop[k] && !closedTest$futilityStop[k]) {
simulatedConditionalPower[k + 1, i] <- simulatedConditionalPower[k + 1, i] +
stageResults$conditionalPowerPerStage[k]
}
}
iterations[k, i] <- iterations[k, i] + 1
for (p in 1:2^(gMax - 1)) {
if (!is.na(stageResults$subjectsPerStage[p, k])) {
simulatedSubjectsPerStage[k, i, p] <- simulatedSubjectsPerStage[k, i, p] +
stageResults$subjectsPerStage[p, k]
}
}
for (g in 1:gMax) {
dataIterationNumber[index] <- j
dataStageNumber[index] <- k
dataPopulationNumber[index] <- g
dataEffect[index] <- i
dataSubjectsPopulation[index] <- stageResults$populationSubjectsPerStage[g, k]
dataNumberOfSubjects[index] <- round(sum(stageResults$subjectsPerStage[, k], na.rm = TRUE), 1)
dataNumberOfCumulatedSubjects[index] <- sum(stageResults$subjectsPerStage[, 1:k], na.rm = TRUE)
dataRejectPerStage[index] <- closedTest$rejected[g, k]
dataTestStatistics[index] <- stageResults$testStatistics[g, k]
dataSuccessStop[index] <- closedTest$successStop[k]
if (k < kMax) {
dataFutilityStop[index] <- closedTest$futilityStop[k]
dataConditionalCriticalValue[index] <- stageResults$conditionalCriticalValue[k]
dataConditionalPowerAchieved[index + 1] <- stageResults$conditionalPowerPerStage[k]
}
dataEffectEstimate[index] <- stageResults$overallEffectSizes[g, k]
dataPValuesSeparate[index] <- closedTest$separatePValues[g, k]
index <- index + 1
}
if (!rejectAtSomeStage && any(closedTest$rejected[, k] &
closedTest$selectedPopulations[1:gMax, k] | rejectedPopulationsBefore)) {
simulatedRejectAtLeastOne[i] <- simulatedRejectAtLeastOne[i] + 1
rejectAtSomeStage <- TRUE
}
if ((k < kMax) && (closedTest$successStop[k] || closedTest$futilityStop[k])) {
# rejected hypotheses remain rejected also in case of early stopping
simulatedRejections[(k + 1):kMax, i, ] <- simulatedRejections[(k + 1):kMax, i, ] +
matrix(
(closedTest$rejected[, k] &
closedTest$selectedPopulations[1:gMax, k] | rejectedPopulationsBefore),
kMax - k, gMax,
byrow = TRUE
)
break
}
rejectedPopulationsBefore <- closedTest$rejected[, k] &
closedTest$selectedPopulations[1:gMax, k] | rejectedPopulationsBefore
}
}
simulatedSubjectsPerStage[is.na(simulatedSubjectsPerStage)] <- 0
simulatedSubjectsPerStage[, i, ] <- simulatedSubjectsPerStage[, i, ] / iterations[, i]
if (kMax > 1) {
simulatedRejections[2:kMax, i, ] <- simulatedRejections[2:kMax, i, ] - simulatedRejections[1:(kMax - 1), i, ]
stopping <- cumsum(simulatedSuccessStopping[1:(kMax - 1), i] +
simulatedFutilityStopping[, i]) / maxNumberOfIterations
expectedNumberOfSubjects[i] <- sum(simulatedSubjectsPerStage[1, i, ] + t(1 - stopping) %*%
simulatedSubjectsPerStage[2:kMax, i, ])
} else {
expectedNumberOfSubjects[i] <- sum(simulatedSubjectsPerStage[1, i, ])
}
}
simulatedConditionalPower[1, ] <- NA_real_
if (kMax > 1) {
simulatedConditionalPower[2:kMax, ] <- as.matrix(simulatedConditionalPower[2:kMax, ] / iterations[2:kMax, ])
}
simulationResults$rejectAtLeastOne <- simulatedRejectAtLeastOne / maxNumberOfIterations
simulationResults$numberOfPopulations <- simulatedNumberOfPopulations / iterations
simulationResults$selectedPopulations <- simulatedSelections / maxNumberOfIterations
simulationResults$rejectedPopulationsPerStage <- simulatedRejections / maxNumberOfIterations
simulationResults$successPerStage <- simulatedSuccessStopping / maxNumberOfIterations
simulationResults$futilityPerStage <- simulatedFutilityStopping / maxNumberOfIterations
simulationResults$futilityStop <- base::colSums(simulatedFutilityStopping / maxNumberOfIterations)
if (kMax > 1) {
simulationResults$earlyStop <- simulationResults$futilityPerStage +
simulationResults$successPerStage[1:(kMax - 1), ]
simulationResults$conditionalPowerAchieved <- simulatedConditionalPower
}
simulationResults$sampleSizes <- simulatedSubjectsPerStage
simulationResults$expectedNumberOfSubjects <- expectedNumberOfSubjects
simulationResults$iterations <- iterations
if (!all(is.na(simulationResults$conditionalPowerAchieved))) {
simulationResults$.setParameterType("conditionalPowerAchieved", C_PARAM_GENERATED)
}
if (any(simulationResults$rejectedPopulationsPerStage < 0)) {
stop(C_EXCEPTION_TYPE_RUNTIME_ISSUE, "internal error, simulation not possible due to numerical overflow")
}
data <- data.frame(
iterationNumber = dataIterationNumber,
stageNumber = dataStageNumber,
populationNumber = dataPopulationNumber,
effect = dataEffect,
numberOfSubjects = dataNumberOfSubjects,
numberOfCumulatedSubjects = dataNumberOfCumulatedSubjects,
subjectsPopulation = dataSubjectsPopulation,
effectEstimate = dataEffectEstimate,
testStatistics = dataTestStatistics,
pValue = dataPValuesSeparate,
conditionalCriticalValue = round(dataConditionalCriticalValue, 6),
conditionalPowerAchieved = round(dataConditionalPowerAchieved, 6),
rejectPerStage = dataRejectPerStage,
successStop = dataSuccessStop,
futilityPerStage = dataFutilityStop
)
data <- data[!is.na(data$effectEstimate), ]
simulationResults$.data <- data
return(simulationResults)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.