Nothing
## Exact simulation for the 3+3 design
setOldClass(c("exact","three_plus_three_selector_factory","tox_selector_factory","selector_factory"))
#' A wrapper function supporting exact simulation of dose-escalation trials.
#'
#' Implemented currently only for the (most?) common variant of the 3+3 design,
#' which requires that at least 6 patients be treated at a dose before it may be
#' declared as \sQuote{the} MTD.
#'
#' @param selector_factory An object of type
#' \code{\link[escalation:get_three_plus_three]{three_plus_three_selector_factory}},
#' with `allow_deescalation = TRUE`.
#'
#' @details
#' In any given realization of a 3+3 design, each of the \eqn{D} prespecified doses
#' will enroll 0, 1 or 2 cohorts, each with 3 patients. Each cohort will result in
#' a tally of 0--3 dose-limiting toxicities (DLTs), and these may be recorded in a
#' \eqn{2 \times D}{2 x D} matrix. Moreover, the 3+3 dose-escalation rules allow for
#' only one path through any such matrix. For example, the matrix
#' \preformatted{
#' d
#' c D1 D2 D3 D4
#' 1 0 1 2 NA
#' 2 NA 0 NA NA
#' }
#' represents the path in a 4-dose 3+3 trial, where the following events occur:
#' \enumerate{
#' \item Initial cohort at \eqn{d=1} results 0/3
#' \item Escalation to \eqn{d=2} results 1/3
#' \item Additional cohort at \eqn{d=2} results 0/3 for net 1/6 at this dose
#' \item Escalation to \eqn{d=3} results 2/3; MTD declared as \eqn{d=1}.
#' }
#' (Indeed, as you may verify at the R prompt, the above matrix is the 262nd of 442
#' such paths enumerated comprehensively in the \eqn{2 \times 4 \times 442}{2 x 4 x 442}
#' array `precautionary:::T[[4]]`.)
#'
#' As detailed in Norris 2020c (below), these matrices may be used to construct simple
#' matrix computations that altogether eliminate the need for discrete-event simulation
#' of the 3+3 design. For each \eqn{D = 2,...,8}, the `precautionary` package has
#' pretabulated a \eqn{J \times 2D}{J x 2D} matrix `precautionary:::U[[D]]` and
#' \eqn{J}-vector `precautionary:::b[[D]]` such that the \eqn{J}-vector \eqn{\pi}
#' of path probabilities is given by:
#' \deqn{
#' log(\pi) = b + U [log(p), log(q)]',
#' }
#' where \eqn{p} is the \eqn{D}-vector of DLT probabilities at the prespecified
#' doses, and \eqn{q \equiv 1-p}{q := 1-p} is its complement. See Eq. (4) of
#' Norris (2020c).
#'
#' For details on the enumeration itself, please see the Prolog code in directory
#' `exec/` of the installed package.
#' @references
#' Norris DC. What Were They Thinking? Pharmacologic priors implicit in a choice
#' of 3+3 dose-escalation design. arXiv:2012.05301 \[stat.ME\]. December 2020.
#' <https://arxiv.org/abs/2012.05301>
#' @examples
#' # Run an exact version of the simulation from FDA-proactive vignette
#' design <- get_three_plus_three(
#' num_doses = 6
#' , allow_deescalate = TRUE)
#' old <- options(
#' dose_levels = c(2, 6, 20, 60, 180, 400)
#' , ordinalizer = function(MTDi, r0 = 1.5)
#' MTDi * r0 ^ c(Gr1=-2, Gr2=-1, Gr3=0, Gr4=1, Gr5=2)
#' )
#' mtdi_gen <- hyper_mtdi_lognormal(CV = 0.5
#' ,median_mtd = 180
#' ,median_sdlog = 0.6
#' ,units="ng/kg/week")
#' exact(design) %>% simulate_trials(
#' num_sims = 1000
#' , true_prob_tox = mtdi_gen) -> EXACT
#' summary(EXACT)$safety
#' if (interactive()) { # runs too long for CRAN servers
#' # Compare with discrete-event-simulation trials
#' design %>% simulate_trials(
#' num_sims = 1000
#' , true_prob_tox = mtdi_gen) -> DISCRETE
#' summary(DISCRETE)$safety[,]
#' # Note the larger MCSEs in this latter simulation, reflecting combined noise
#' # from hyperprior sampling and the nested discrete-event trial simulations.
#' # The MCSE of the former simulation is purely from the hyperprior sampling,
#' # since the nested trial simulation is carried out by an exact computation.
#' }
#' options(old)
#' @export
exact <- function(selector_factory) {
stopifnot("Class 'exact' applies only to a three_plus_three_selector_factory"
= is(selector_factory,"three_plus_three_selector_factory"))
stopifnot("Exact sim is currently implemented only for 3+3 with 'allow_deescalate=TRUE'"
= selector_factory$allow_deescalate)
prependClass("exact", selector_factory)
}
# Calculate the G matrix from an mtdi_dist
# NB: The ... may be used to pass r0 thru to ordinalizer
G <- function(mtdi_dist, ordinalizer = getOption("ordinalizer"), ...) {
# 1. Get a (D x 5) matrix of doses rescaled by r0^(3-g), g in 1:5
# 2. Get the CDF values on this matrix
# Note that the following assumes ordinalizers multiply by a 5-vector
# of the form c(r_2, r_1, 1, r1, r2), where r_2 < r_1 < 1 < r1 < r2.
r0_spread <- 1/ordinalizer(1, ...)
G5 <- mtdi_dist@dist$cdf(getOption("dose_levels") %o% r0_spread)
# 3. Attach trivial columns of 1's and 0's; label for clarity
G7 <- cbind(None=rep(1,nrow(G5)), G5, rep(0,nrow(G5)))
dimnames(G7)[1] <- list(D = paste(getOption('dose_levels'), mtdi_dist@units))
# The (now) 7 columns give the probabilities of toxicity Gr=0+, Gr=1+, ..., Gr=6+.
# The trivial leftmost and rightmost columns assist with difference-taking below.
# 4. By taking differences, we obtain probabilities for the grades *categorically*
G <- G7[,1:6] - G7[,2:7]
# 5. Form the (2D x 6) blocked matrix (note, this remains a very small matrix!)
Z <- matrix(0, nrow = nrow(G), ncol = 3)
dimnames(Z) <- dimnames(G[,1:3])
G <- rbind(cbind(Z, G[,4:6]),
cbind(G[,1:3], Z))
# 6. Normalize the rows
G <- G / rowSums(G)
# If we do not bound CV away from zero, then it is possible to draw an mtdi_dist
# that e.g. makes non-DLTs *impossible* at some high dose -- i.e., all toxicities
# are of grade 3 or higher. In such extreme cases, the normalizing factor rowSums(G)
# will have zero in components corresponding to such doses. (The converse may also
# occur, in which low doses *cannot* cause a DLT.) This results in NaN's in the G
# matrix, which cascade into results despite multiplication by zero probabilities
# (corresponding to impossible events) in the U matrix. To avoid this, we set these
# elements of the matrix to 0.
G[is.nan(G)] <- 0
return(G)
}
#' @importFrom escalation prob_recommend
prob_recommend.exact <- function(x, ...) {
prob_recs <- NextMethod()
names(prob_recs)[-1] <- paste0(x$dose_levels, x$dose_units)
prob_recs
}
#' @importFrom escalation prob_administer
prob_administer.exact <- function(x, ...) {
prob_admin <- NextMethod()
names(prob_admin) <- paste0(x$dose_levels, x$dose_units)
prob_admin
}
#' Specialize print method for objects of class \code{\link[escalation]{simulations}}
#'
#' @param x An object of class c("exact","precautionary","simulations")
#'
#' @param ... Additional arguments; ignored
#'
#' @importFrom escalation num_patients num_tox trial_duration
print.exact <- function(x, ...) {
cat('Number of iterations:', length(x$fits), '\n')
cat('\n')
cat('Number of doses:', length(x$dose_levels), '\n')
cat('\n')
cat('True probability of toxicity:\n')
print(x$true_prob_tox, digits = 3)
cat('\n')
cat('Probability of recommendation:\n')
print(prob_recommend(x), digits = 3)
cat('\n')
cat('Probability of administration:\n')
print(prob_administer(x), digits = 3)
cat('\n')
cat('Sample size:\n')
print(summary(num_patients(x)))
cat('\n')
cat('Total toxicities:\n')
print(summary(num_tox(x)))
cat('\n')
cat('Trial duration:\n')
print(summary(trial_duration(x)))
cat('\n')
}
#' Specialize print method defined for class \code{\link[escalation]{simulations}}
#'
#' @param x An object of class c("hyper","precautionary","simulations")
#'
#' @param ... Additional arguments; ignored
#'
#' @importFrom escalation num_patients num_tox trial_duration
print.hyper <- function(x, ...) {
cat('Number of iterations:', length(x$fits), '\n')
cat('\n')
cat('Number of doses:', length(x$dose_levels), '\n')
cat('\n')
cat('Average probability of toxicity:\n')
print(x$avg_prob_tox, digits = 3)
cat('\n')
cat('Probability of recommendation:\n')
print(prob_recommend(x), digits = 3)
cat('\n')
cat('Probability of administration:\n')
print(prob_administer(x), digits = 3)
cat('\n')
cat('Sample size:\n')
print(summary(num_patients(x)))
cat('\n')
cat('Total toxicities:\n')
print(summary(num_tox(x)))
cat('\n')
cat('Trial duration:\n')
print(summary(trial_duration(x)))
cat('\n')
}
#' Convert an object of class c('exact','precautionary','simulations') to a data.table
#'
#' TODO: Actually implement this!
#' @param x An object of class c('precautionary','simulations')
#'
#' @param keep.rownames Unused; retained for S3 generic/method consistency
#' @param ordinalizer If not NULL, this is a function mapping the threshold
#' dose ('MTDi') at which an individual experiences a binary toxicity (as
#' recognized by the dose-escalation design) to a named vector giving dose
#' thresholds for multiple grades of toxicity. The names of this vector will
#' be taken as designations of the toxicity grades.
#' @param ... Additional parameters passed to the `ordinalizer`
#'
#' @export
as.data.table.exact <- function(x, keep.rownames = FALSE
, ordinalizer = getOption('ordinalizer')
, ...) {
extractor <- ifelse(is(x$fits[[1]][[1]]$fit, "derived_dose_selector")
,function(.) .[[1]]$fit$parent$outcomes
,function(.) .[[1]]$fit$outcomes
)
ensemble <- rbindlist(lapply(x$fits, extractor), idcol = "rep")
if( is.null(ordinalizer) )
return(ensemble)
# TODO: Do add these columns to 'ensemble', so that the whole table
# may later be inspected by user to improve understanding.
MTDig <- t(sapply(ensemble$MTDi, ordinalizer, ...))
if( is.null(colnames(MTDig)) ) {
warning("Ordinalizer returns unnamed vector; using default names for toxicity grades.")
colnames(MTDig) <- paste("Grade", 1:ncol(MTDig))
}
tox_grades <- colnames(MTDig)
# Compare with actual dose to obtain toxicity grade indicator matrix
tox_ind <- ( x$dose_levels[ensemble$dose] > MTDig )
# Tally the thresholds crossed to obtain integer toxgrade
ensemble$toxgrade <- rowSums(tox_ind)
# Convert toxgrade to an ordered factor Tox
ensemble$Tox <- ordered(ensemble$toxgrade+1
, levels=seq(1+length(tox_grades))
, labels=c('None', tox_grades))
ensemble
}
#' Summarize an exact treatment of a dose-escalation design
#'
#' Algorithmic (or 'rule-based') dose-escalation designs admit exact computation
#' of their outcomes. This method summarizes such an exact treatment in a manner
#' roughly parallel to that of `summary.precautionary`.
#'
#' @param object An object of class 'exact'
#'
#' @param ordinalizer An ordinalizer function
#' @param ... Additional parameters passed to the ordinalizer
#'
#' @importFrom dplyr mutate rename_with select everything
#' @importFrom stats xtabs addmargins sd var
#' @importFrom rlang .data
#' @export
summary.exact <- function(object, ordinalizer = getOption('ordinalizer'), ...) {
# NB: We are severing the call chain, with a special implementation here
##summary <- NextMethod()
# TODO: Consider implementing the 'escalation' & 'toxTab' components.
# (I believe both are obtainable in the exact formulation.)
# I need to deal with the 'hyper' case separately.
summary <- list(escalation = "not yet implemented",
safety = NULL, # to be calculated ...
toxTab = "not yet implemented")
if (is(object, 'hyper')) {
toxTab <- object$hyper$safety %>%
addmargins(margin = 2, FUN = list(Total=sum))
summary$safety <- rbind(
"Expected participants" = colMeans(toxTab)
,"MCSE" = apply(toxTab, MARGIN = 2, FUN = sd) / sqrt(nrow(toxTab))
)
} else {
summary$safety <- object$safety %>%
addmargins(margin = 2, FUN = list(Total=sum))
}
summary
}
# TODO: Incorporate this logic to map 3+3 sim outcomes indices j in T[[D]][,,j]
#' @importFrom stats aggregate
haystack <- function(sims) {
fit <- sims$fits[[1]][[1]]$fit
m <- function(fit) {
oc <- fit$outcomes
ag <- aggregate(oc$tox, by=oc[,c("cohort","dose")], FUN=sum)
ag$cohort <- NULL # drop this column
# For any doses appearing just once, add an NA dose
jo <- which(tabulate(ag$dose) == 1) # jo = 'just once'
if (length(jo))
ag <- rbind(ag, data.frame(dose = jo, x = NA))
# For doses that never appeared, add 2 NAs
doses <- 1:fit$num_doses
nd <- doses[doses > max(ag$dose)]
if (length(nd))
ag <- rbind(ag, data.frame(dose = rep(nd, 2), x = NA))
m <- matrix(ag[order(ag$dose),"x"], nrow = 2)
dimnames(m) <- list(c = 1:2, d = paste0("D",doses))
m
}
j <- function(m) {
D <- ncol(m)
J <- dim(T[[D]])[3]
which(sapply(1:J, function(j)
identical(m, T[[D]][,,j])))
}
j(m(fit))
}
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.