Nothing
# Try helper functions -----
#' Tries to Evaluate an Expression and Fails with \code{NA}
#'
#' The evaluation fails with \code{NA} by default, but it is also able to fail with other values.
#'
#' @param expr Expression to be evaluated.
#' @param value Return value if there is an error, default is \code{NA_real_}.
#'
#' @return Returns the evaluation of the expression, or \code{value} if it doesn't work out.
tryOrFailWithNA <- function(expr, value=NA_real_) {
tryCatch(
error=function(cnd) value,
expr
)
}
#' Checks Whether a Vector of Object Inherits from the Class 'try-error'
#'
#' Checks whether any of the provided objects contains a try error.
#'
#' @param ... objects that need testing.
#'
#' @return Returns \code{TRUE} if there's some object that's a try-error, \code{FALSE} when all objects are
#' not try-errors.
#'
#' @export
#'
#' @examples
#' x <- 1
#' y <- "a"
#' z <- try(integrate(exp, -Inf, Inf))
#' isTryError(x, y)
#' isTryError(x, y, z)
isTryError <- function(...) {
obj <- list(...)
tryErrorFunc <- function(x){inherits(x, "try-error")}
result <- purrr::some(obj, .p=tryErrorFunc)
return(result)
}
#' Helper function: Get all arguments as entered by the user
#'
#' @return a list of variable names of class "call" that can be changed into names
getArgs <- function() {
as.list(match.call(definition = sys.function(-1),
call = sys.call(-1)))[-1]
}
#' Helper function: Get all names as entered by the user
#'
#' @param list list from which the element needs retrieving
#' @param name character string, name of the item that need retrieving
#'
#' @return returns a character string
extractNameFromArgs <- function(list, name) {
result <- list[[name]]
if (inherits(result, "call"))
result <- as.character(as.expression(result))
return(result)
}
# Check Consistency function --------
#' Checks consistency between the sided of the hypothesis and the minimal clinically relevant effect size
#' or safe test defining parameter. Throws an error if the one-sided hypothesis is incongruent with the
#'
#' @inheritParams designSafeZ
#' @param paramToCheck numeric. Either a named safe test defining parameter such as phiS, or thetaS, or a
#' minimal clinically relevant effect size called with a non-null esMinName name
#' @param esMinName provides the name of the effect size. Either "meanDiffMin" for the z-test, "deltaMin" for
#' the t-test, or "hrMin" for the logrank test
#' @param paramDomain Domain of the paramToCheck, typically, positiveNumbers. Default \code{NULL}
#'
#' @return paramToCheck after checking, perhaps with a change in sign
checkAndReturnsEsMinParameterSide <- function(paramToCheck, alternative=c("twoSided", "greater", "less"),
esMinName=c("noName", "meanDiffMin", "phiS",
"deltaMin", "deltaS",
"hrMin", "thetaS", "deltaTrue"),
paramDomain=NULL) {
# TODO(Alexander): Remove in v0.9.0
#
if (length(alternative)==1 && alternative=="two.sided") {
warning('The option alternative="two.sided" is deprecated;',
'Please use alternative="twoSided" instead')
alternative <- "twoSided"
}
alternative <- match.arg(alternative)
paramDomain <- match.arg(paramDomain)
esMinName <- match.arg(esMinName)
if (alternative == "twoSided") {
if (esMinName %in% c("meanDiffMin", "deltaMin", "deltaTrue"))
return(abs(paramToCheck))
return(paramToCheck)
}
if (esMinName=="noName")
paramName <- NULL
else
paramName <- esMinName
if (is.null(paramName)) {
paramName <- "the safe test defining parameter"
hypParamName <- "test relevant parameter"
paramDomain <- "unknown"
} else if (paramName=="phiS" || esMinName=="meanDiffMin") {
hypParamName <- "meanDiff"
paramDomain <- "realNumbers"
} else if (paramName=="deltaS" || esMinName=="deltaMin" || esMinName=="deltaTrue") {
hypParamName <- "delta"
paramDomain <- "realNumbers"
} else if (paramName=="thetaS" || esMinName=="hrMin") {
hypParamName <- "theta"
paramDomain <- "positiveNumbers"
} else {
hypParamName <- "testRelevantParameter"
}
if (paramDomain=="unknown") {
nullValue <- "nullValue"
if (alternative=="greater" && paramToCheck < 0) {
warning('The safe test defining parameter is incongruent with alternative "greater". ',
"This safe test parameter is made positive to compare H+: ",
"test-relevant parameter > 0 against H0 : test-relevant parameter = 0")
paramToCheck <- -paramToCheck
}
if (alternative=="less" && paramToCheck > 0) {
warning('The safe test defining parameter is incongruent with alternative "less". ',
"This safe test parameter is made positive to compare H-: ",
"test-relevant parameter < 0 against H0 : test-relevant parameter = 0")
paramToCheck <- -paramToCheck
}
} else if (paramDomain=="realNumbers") {
nullValue <- 0
if (alternative=="greater" && paramToCheck < 0) {
warning(paramName, ' incongruent with alternative "greater". ',
paramName, " set to -", paramName, " > 0 in order to compare H+: ",
hypParamName, " > 0 against H0 : ", hypParamName, " = 0")
paramToCheck <- -paramToCheck
}
if (alternative=="less" && paramToCheck > 0) {
warning(paramName, ' incongruent with alternative "greater". ',
paramName, " set to -", paramName, " < 0 in order to compare H-: ",
hypParamName, " < 0 against H0 : ", hypParamName, " = 0")
paramToCheck <- -paramToCheck
}
} else if (paramDomain=="positiveNumbers") {
if (alternative=="greater" && paramToCheck < 1) {
warning(paramName, ' incongruent with alternative "greater". ',
paramName, " set to 1/", paramName, " > 1 in order to compare H+: ",
hypParamName, " > 1 against H0 : ", hypParamName, " = 1")
paramToCheck <- 1/paramToCheck
}
if (alternative=="less" && paramToCheck > 1) {
warning(paramName, ' incongruent with alternative "greater". ',
paramName, " set to 1/", paramName, " < 1 in order to compare H-: ",
hypParamName, " < 1 against H0 : ", hypParamName, " = 1")
paramToCheck <- 1/paramToCheck
}
}
return(paramToCheck)
}
#' Check consistency between nPlan and the testType for one and two-sample z and t-tests
#'
#' @inheritParams designSafeZ
#'
#' @return nPlan a vector of sample sizes of length 1 or 2
#'
checkAndReturnsNPlan <- function(nPlan, ratio=1, testType=c("oneSample", "paired", "twoSample")) {
if (testType=="twoSample" && length(nPlan)==1) {
nPlan <- c(nPlan, ratio*nPlan)
warning('testType=="twoSample" specified, but nPlan[2] not provided. nPlan[2] = ratio*nPlan[1], that is, ',
nPlan[2], '.')
} else if (testType=="paired" && length(nPlan)==1) {
nPlan <- c(nPlan, nPlan)
warning('testType=="paired" specified, but nPlan[2] not provided. nPlan[2] set to nPlan[1].')
} else if (testType=="oneSample" && length(nPlan)==2) {
nPlan <- nPlan[1]
warning('testType=="oneSample" specified, but two nPlan[2] provided, which is ignored.')
}
return(nPlan)
}
# Plot helper -----
#' Sets 'safestats' Plot Options and Returns the Current Plot Options.
#'
#' @param ... further arguments to be passed to or from methods.
#'
#' @return Returns a list with the user specified plot options.
#'
#' @export
#'
#' @examples
#' oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes()
#' graphics::plot(1:10, 1:10)
#' setPar <- graphics::par(oldPar)
setSafeStatsPlotOptionsAndReturnOldOnes <- function(...) {
oldPar <- graphics::par(no.readonly = TRUE)
graphics::par(cex.main=1.5, mar=c(5, 6, 4, 4)+0.1, mgp=c(3.5, 1, 0), cex.lab=1.5,
font.lab=2, cex.axis=1.3, bty="n", las=1, ...)
return(oldPar)
}
# Vignette helpers ---------
#' Plots the Histogram of Stopping Times
#'
#' Helper function to display the histogram of stopping times.
#'
#' @param safeSim A safeSim object, returned from \code{\link{replicateTTests}}.
#' @param nPlan numeric > 0, the planned sample size(s).
#' @param deltaTrue numeric, that represents the true underlying standardised effect size delta.
#' @param showOnlyNRejected logical, when \code{TRUE} discards the cases that did not reject.
#' @param nBin numeric > 0, the minimum number of bins in the histogram.
#' @param ... further arguments to be passed to or from methods.
#'
#' @return a histogram object, and called for its side-effect to plot the histogram.
#'
#' @export
#'
#' @examples
#'# Design safe test
#' alpha <- 0.05
#' beta <- 0.20
#' designObj <- designSafeT(1, alpha=alpha, beta=beta)
#'
#' # Design frequentist test
#' freqObj <- designFreqT(1, alpha=alpha, beta=beta)
#'
#' # Simulate under the alternative with deltaTrue=deltaMin
#' simResults <- replicateTTests(nPlan=designObj$nPlan, deltaTrue=1, parameter=designObj$parameter,
#' nPlanFreq=freqObj$nPlan)
#'
#' plotHistogramDistributionStoppingTimes(
#' simResults$safeSim, nPlan = simResults$nPlan,
#' deltaTrue = simResults$deltaTrue)
plotHistogramDistributionStoppingTimes <- function(safeSim, nPlan, deltaTrue, showOnlyNRejected=FALSE, nBin=25L, ...) {
if(showOnlyNRejected) {
dataToPlot <- safeSim[["allRejectedN"]]
} else {
dataToPlot <- safeSim[["allN"]]
}
nStep <- floor(nPlan/nBin)
if (nStep == 0)
nStep <- 1
maxLength <- ceiling(nPlan/nStep)
mainTitle <- bquote(~"Spread of stopping times when true delta " == .(round(deltaTrue,2)))
oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes()
on.exit(graphics::par(oldPar))
graphics::hist(dataToPlot,
breaks = nStep*seq.int(maxLength),
xlim = c(0, max(safeSim[["allN"]])),
xlab = "stopping time (n collected)",
main = mainTitle,
col = "lightgrey", ...)
}
#' Selectively Continue Experiments that Did Not Lead to a Null Rejection for a (Safe) T-Test
#'
#' Helper function used in the vignette.
#'
#' @inheritParams replicateTTests
#' @param oldValues vector of e-values or p-values.
#' @param valuesType character string either "eValues" or "pValues".
#' @param designObj a safeDesign object obtained from \code{\link{designSafeT}}, or \code{NULL}
#' if valuesType equal "pValues".
#' @param oldData a list of matrices with names "dataGroup1" and "dataGroup2".
#' @param n1Extra integer, that defines the additional number of samples of the first group. If
#' \code{NULL} and valuesType equals "eValues", then n1Extra equals \code{designObj$nPlan[1]}.
#' @param n2Extra optional integer, that defines the additional number of samples of the second group.
#' If \code{NULL}, and valuesType equals "eValues", then n2Extra equals \code{designObj$nPlan[2]}.
#' @param moreMainText character, additional remarks in the title of the histogram.
#'
#' @return a list that includes the continued s or p-values based on the combined data, and a list of the
#' combined data.
#' @export
#'
#' @examples
#' alpha <- 0.05
#' mIter <- 1000L
#'
#' designObj <- designSafeT(deltaMin=1, alpha=alpha, beta=0.2, nSim=100)
#' oldData <- generateNormalData(nPlan=designObj$nPlan, deltaTrue=0, nSim=mIter, seed=1)
#'
#' eValues <- vector("numeric", length=mIter)
#'
#' for (i in seq_along(eValues)) {
#' eValues[i] <- safeTTest(x=oldData$dataGroup1[i, ], designObj=designObj)$eValue
#' }
#'
#' # First run: 8 false null rejections
#' sum(eValues > 1/alpha)
#'
#' continuedSafe <- selectivelyContinueTTestCombineData(
#' oldValues=eValues, designObj=designObj, oldData=oldData,
#' deltaTrue=0, seed=2)
#'
#' # Second run: 1 false null rejections
#' sum(continuedSafe$newValues > 1/alpha)
#'
#' # Third run: 0 false null rejections
#' eValues <- continuedSafe$newValues
#' oldData <- continuedSafe$combinedData
#' continuedSafe <- selectivelyContinueTTestCombineData(
#' oldValues=eValues, designObj=designObj, oldData=oldData,
#' deltaTrue=0, seed=3)
#' sum(continuedSafe$newValues > 1/alpha)
selectivelyContinueTTestCombineData <- function(oldValues, valuesType=c("eValues", "pValues"), designObj=NULL,
alternative=c("twoSided", "greater", "less"),
oldData, deltaTrue, alpha=NULL,
n1Extra=NULL, n2Extra=NULL, seed=NULL, paired=FALSE,
muGlobal=0, sigmaTrue=1, moreMainText="") {
valuesType <- match.arg(valuesType)
# TODO(Alexander): Remove in v0.9.0
#
if (length(alternative)==1 && alternative=="two.sided") {
warning('The option alternative="two.sided" is deprecated;',
'Please use alternative="twoSided" instead')
alternative <- "twoSided"
}
alternative <- match.arg(alternative)
if (valuesType=="pValues") {
if (is.null(alpha)) {
stop('For valuesType="pValues" an alpha is required')
}
if (is.null(n1Extra) && is.null(n2Extra)) {
stop("Can't sample extra data without being specified the number of extra samples")
}
} else if (valuesType=="eValues") {
if (is.null(designObj)) {
stop("Can't continue safe test analysis without a design object")
}
alpha <- designObj[["alpha"]]
if (is.null(n1Extra) && is.null(n2Extra)) {
n1Extra <- designObj[["nPlan"]][1]
n2Extra <- designObj[["nPlan"]][2]
}
}
if (is.null(oldData[["dataGroup1"]]))
stop("Can't combine data from old experiment without a matrix of old data referred to as oldData$dataGroup1")
result <- list("oldValues"=oldValues, "newValues"=NULL, "valuesType"=valuesType, "combinedData"=NULL, "deltaTrue"=deltaTrue,
"cal"=sys.call())
if (valuesType=="eValues") {
notRejectedIndex <- which(oldValues <= 1/alpha)
} else if (valuesType=="pValues") {
notRejectedIndex <- which(oldValues >= alpha)
}
if (length(notRejectedIndex)==0)
stop("All experiments led to a null rejection. Nothing to selectively continue.")
oldDataGroup1 <- oldData[["dataGroup1"]][notRejectedIndex, ]
oldDataGroup2 <- oldData[["dataGroup2"]][notRejectedIndex, ]
if (is.null(n2Extra) || is.na(n2Extra))
tempNPlan <- n1Extra
else
tempNPlan <- c(n1Extra, n2Extra)
newData <- generateNormalData("nPlan"=tempNPlan,
"deltaTrue"=deltaTrue, "nSim"=length(notRejectedIndex),
"paired"=paired, "seed"=seed,
"muGlobal"=muGlobal, "sigmaTrue"=sigmaTrue)
dataGroup1 <- cbind(oldDataGroup1, newData[["dataGroup1"]])
dataGroup2 <- cbind(oldDataGroup2, newData[["dataGroup2"]])
newValues <- vector("numeric", length(notRejectedIndex))
if (valuesType=="eValues") {
for (i in seq_along(newValues)) {
newValues[i] <- try(safeTTest("x"=dataGroup1[i, ], "y"=dataGroup2[i, ],
"designObj"=designObj, "alternative"=alternative,
"paired"=paired)$eValue)
}
minX <- log(min(oldValues, newValues))
maxX <- log(max(oldValues, newValues))
yOld <- log(oldValues[notRejectedIndex])
yNew <- log(newValues)
xLabText <- "log(eValues)"
mainText <- paste("Histogram of e-values", moreMainText)
threshValue <- log(1/alpha)
} else if (valuesType=="pValues") {
for (i in seq_along(newValues)) {
alternativeFreq <- alternative
if (alternative=="twoSided")
alternativeFreq <- "two.sided"
newValues[i] <- t.test("x"=dataGroup1[i, ], "y"=dataGroup2[i, ],
"alternative"=alternativeFreq, "paired"=paired)$p.value
}
minX <- 0
maxX <- 1
yOld <- oldValues[notRejectedIndex]
yNew <- newValues
xLabText <- "p-values"
mainText <- "Histogram of p-values"
threshValue <- alpha
}
oldHist <- graphics::hist("x"=yOld, plot=FALSE)
newHist <- graphics::hist("x"=yNew, plot=FALSE)
yMax <- max(oldHist[["counts"]], newHist[["counts"]])
oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes()
on.exit(graphics::par(oldPar))
graphics::plot(oldHist, "xlim"=c(minX, maxX), ylim=c(0, yMax), "col"="blue",
"density"=20, "angle"=45, "xlab"=xLabText, "main"=mainText)
graphics::plot(newHist, "add"=TRUE, "col"="red", "density"=20, "angle"=-45)
graphics::abline("v"=threshValue, "col"="grey", "lwd"=2, "lty"=2)
result[["newValues"]] <- newValues
result[["combinedData"]] <-list("dataGroup1"=dataGroup1, "dataGroup2"=dataGroup2)
return(result)
}
#' Generate Survival Data which Can Be Analysed With the `survival` Package
#'
#'
#' @param nP integer > 0 representing the number of of patients in the placebo group.
#' @param nT integer > 0 representing the number of of patients in the treatment group.
#' @param alpha numeric > 0, representing the shape parameter of the Weibull distribution.
#' If alpha=1, then data are generated from the exponential, i.e., constant hazard. For alpha > 1
#' the hazard increases, if alpha < 1, the hazard decreases.
#' @param lambdaP The (relative) hazard of the placebo group.
#' @param lambdaT The (relative) hazard of the treatment group.
#' @param seed A seed number.
#' @param nDigits numeric, the number of digits to round of the random time to
#' @param startTime numeric, adds this to the random times. Default 1, so the startTime is not 0, which
#' is the start time of \code{\link[stats]{rweibull}}.
#' @param endTime The endtime of the experiment.
#' @param orderTime logical, if \code{TRUE} then put the data set in increasing order
#' @param competeRatio The ratio of the data that is due to competing risk.
#'
#' @return A data set with time, status and group.
#' @export
#'
#' @examples
#' generateSurvData(800, 800, alpha=1, lambdaP=0.008, lambdaT=0.008/2)
generateSurvData <- function(nP, nT, alpha=1, lambdaP, lambdaT, seed=NULL, nDigits=0,
startTime=1, endTime=180, orderTime=TRUE, competeRatio=0) {
stopifnot(competeRatio >=0, competeRatio < 1, is.numeric(startTime), is.numeric(endTime))
set.seed(seed)
data <- list()
if (competeRatio==0) {
data[["time"]] <- round(c(stats::rweibull("n" = nP, "shape" = alpha,
"scale" = lambdaP^(-1/alpha)),
stats::rweibull("n" = nT, "shape" = alpha,
"scale" = lambdaT^(-1/alpha))) + startTime,
"digits" = nDigits)
data[["status"]] <- 2 # 2 is death
data[["group"]] <- c(rep("P", times = nP), rep("T", times = nT))
data <- as.data.frame(data)
data[["status"]][data[["time"]] > endTime] <- 1 # 1 is censored
data[["time"]][data[["time"]] > endTime] <- endTime
} else {
moreNP <- ceiling((1+competeRatio)*nP)
moreNT <- ceiling((1+competeRatio)*nT)
data[["time"]] <- round(c(stats::rweibull("n" = moreNP, "shape" = alpha,
"scale" = lambdaP^(-1/alpha)),
stats::rweibull("n" = moreNT, "shape" = alpha,
"scale" = lambdaT^(-1/alpha))) + startTime,
"digits" = nDigits)
data[["status"]] <- 2 # 2 is death
data[["group"]] <- c(rep("P", times = moreNP), rep("T", times = moreNT))
data <- as.data.frame(data)
data[["status"]][data[["time"]] > endTime] <- 0 # 0 is censored
data[["time"]][data[["time"]] > endTime] <- endTime
indexOfDeaths <- which(data[["status"]]==2)
totalNDeaths <- length(indexOfDeaths)
nCompete <- floor(competeRatio*totalNDeaths)
if (nCompete < 1)
nCompete <- 1
indexOfCompete <- sample(indexOfDeaths, nCompete)
data[["status"]][indexOfCompete] <- 1
data[["status"]] <- factor(data[["status"]], 0:2, c("censored", "competing", "death"))
}
if (orderTime)
data <- data[order(data[["time"]]), ]
return(data)
}
#' Helper function to check whether arguments are specified in a function at a higher level and already
#' provided in the design object.
#'
#' @param designObj an object of class "safeDesign".
#' @param ... arguments that need checking.
#'
#' @return Returns nothing only used for its side-effects to produces warnings if needed.
#'
#' @export
#'
#' @examples
#' designObj <- designSafeZ(0.4)
#'
#' checkDoubleArgumentsDesignObject(designObj, "alpha"=NULL, alternative=NULL)
#' # Throws a warning
#' checkDoubleArgumentsDesignObject(designObj, "alpha"=0.4, alternative="d")
checkDoubleArgumentsDesignObject <- function(designObj, ...) {
argsToCheck <- list(...)
for (neem in names(argsToCheck)) {
argument <- argsToCheck[[neem]]
if (!is.null(argument) && argument != designObj[[neem]])
warning("Both a design object and '", neem, "' provided. The '", neem, "' specified by the design ",
"object is used for the test, and the provided '", neem, "' is ignored.")
}
}
# ---------- Boot helpers --------
#' Computes the bootObj for sequential sampling procedures regarding nPlan, beta, the implied target
#'
#' @inheritParams designSafeZ
#' @param values numeric vector. If objType equals "nPlan" or "beta" then values should be stopping times,
#' if objType equals "logImpliedTarget" then values should be eValues.
#' @param nBoot integer > 0 representing the number of bootstrap samples to assess the accuracy of
#' approximation of the power, the planned sample size(s) of the safe test under continuous monitoring.
#' @param nPlan integer > 0 representing the number of planned samples (for the first group).
#' @param objType character string either "nPlan", "nMean", "beta", "betaFromEValues", "expectedStopTime" or "logImpliedTarget".
#'
#' @return bootObj
#' @export
#'
#' @examples
#' computeBootObj(1:100, objType="nPlan", beta=0.3)
computeBootObj <- function(values, beta=NULL, nPlan=NULL, nBoot=1e3L, alpha=NULL,
objType=c("nPlan", "nMean", "beta", "betaFromEValues",
"logImpliedTarget", "expectedStopTime")) {
objType <- match.arg(objType)
if (objType=="beta") {
if (is.null(nPlan) || nPlan <= 0)
stop("Please provide an nPlan > 0")
times <- values
stopifnot(nPlan > 0)
bootObj <- boot::boot(times,
function(x, idx) {
1-mean(x[idx] <= nPlan)
}, R = nBoot)
} else if (objType =="betaFromEValues") {
if (is.null(alpha) || alpha <= 0 || alpha >= 1)
stop("Please provide an alpha in (0, 1)")
eValues <- values
bootObj <- boot::boot(
data = eValues,
statistic = function(x, idx) {
mean(x[idx] >= 1/alpha)
},
R = nBoot
)
} else if (objType=="nPlan") {
if (is.null(beta) || beta <= 0 || beta >= 1)
stop("Please provide a beta in (0, 1)")
times <- values
bootObj <- boot::boot(times,
function(x, idx) {
quantile(x[idx], prob=1-beta, names=FALSE)
}, R = nBoot)
} else if (objType=="nMean") {
if (is.null(nPlan[1]) || nPlan[1] <= 0)
stop("Please provide a positive nPlan")
times <- values
times[times > nPlan[1]] <- nPlan[1]
bootObj <- boot::boot(times,
function(x, idx) {
mean(x[idx])
}, R = nBoot)
} else if (objType=="logImpliedTarget") {
eValues <- values
stopifnot(eValues > 0)
bootObj <- boot::boot(eValues,
function(x, idx) {
mean(log(x[idx]))
}, R = nBoot)
} else if (objType=="expectedStopTime") {
times <- values
bootObj <- boot::boot(times,
function(x, idx) {
mean(x[idx])
}, R = nBoot)
}
bootObj[["bootSe"]] <- sd(bootObj[["t"]])
return(bootObj)
}
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.