#' @include checkmate.R
#' @include Model-methods.R
#' @include Samples-class.R
#' @include Rules-class.R
#' @include helpers.R
#' @include helpers_rules.R
#' @include helpers_broom.R
NULL
# nextBest ----
## generic ----
#' Finding the Next Best Dose
#'
#' @description `r lifecycle::badge("stable")`
#'
#' A function that computes the recommended next best dose based on the
#' corresponding rule `nextBest`, the posterior `samples` from the `model` and
#' the underlying `data`.
#'
#' @param nextBest (`NextBest`)\cr the rule for the next best dose.
#' @param doselimit (`number`)\cr the maximum allowed next dose. If it is an
#' infinity (default), then essentially no dose limit will be applied in the
#' course of dose recommendation calculation.
#' @param samples (`Samples`)\cr posterior samples from `model` parameters given
#' `data`.
#' @param model (`GeneralModel`)\cr model that was used to generate the samples.
#' @param data (`Data`)\cr data that was used to generate the samples.
#' @param ... additional arguments without method dispatch.
#'
#' @return A list with the next best dose recommendation (element named `value`)
#' from the grid defined in `data`, and a plot depicting this recommendation
#' (element named `plot`). In case of multiple plots also an element
#' named `singlePlots` is included. The `singlePlots` is itself a list with
#' single plots. An additional list with elements describing the outcome
#' of the rule can be contained too.
#'
#' @export
#'
setGeneric(
name = "nextBest",
def = function(nextBest, doselimit, samples, model, data, ...) {
if (!missing(doselimit)) {
assert_number(doselimit, lower = 0, finite = FALSE)
}
standardGeneric("nextBest")
},
valueClass = "list"
)
## NextBestMTD ----
#' @describeIn nextBest find the next best dose based on the MTD rule.
#'
#' @aliases nextBest-NextBestMTD
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestMTD.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestMTD",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
# Generate the MTD samples and derive the next best dose.
dose_target_samples <- dose(x = nextBest@target, model, samples, ...)
dose_target <- nextBest@derive(dose_target_samples)
# Round to the next possible grid point.
doses_eligible <- h_next_best_eligible_doses(data@doseGrid, doselimit, data@placebo)
next_dose_level <- which.min(abs(doses_eligible - dose_target))
next_dose <- doses_eligible[next_dose_level]
# Create a plot.
p <- ggplot(
data = data.frame(x = dose_target_samples),
aes(.data$x),
fill = "grey50",
colour = "grey50"
) +
geom_density() +
coord_cartesian(xlim = range(data@doseGrid)) +
geom_vline(xintercept = dose_target, colour = "black", lwd = 1.1) +
geom_text(
data = data.frame(x = dose_target),
aes(.data$x, 0),
label = "Est",
vjust = -0.5,
hjust = 0.5,
colour = "black",
angle = 90
) +
xlab("MTD") +
ylab("Posterior density")
if (is.finite(doselimit)) {
p <- p +
geom_vline(xintercept = doselimit, colour = "red", lwd = 1.1) +
geom_text(
data = data.frame(x = doselimit),
aes(.data$x, 0),
label = "Max",
vjust = -0.5,
hjust = -0.5,
colour = "red",
angle = 90
)
}
p <- p +
geom_vline(xintercept = next_dose, colour = "blue", lwd = 1.1) +
geom_text(
data = data.frame(x = next_dose),
aes(.data$x, 0),
label = "Next",
vjust = -0.5,
hjust = -1.5,
colour = "blue",
angle = 90
)
list(value = next_dose, plot = p)
}
)
## NextBestNCRM ----
#' @describeIn nextBest find the next best dose based on the NCRM method. The
#' additional element `probs` in the output's list contains the target and
#' overdosing probabilities (across all doses in the dose grid) used in the
#' derivation of the next best dose.
#'
#' @aliases nextBest-NextBestNCRM
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestNCRM.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestNCRM",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
# Matrix with samples from the dose-tox curve at the dose grid points.
prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
# Estimates of posterior probabilities that are based on the prob. samples
# which are within overdose/target interval.
prob_overdose <- colMeans(h_in_range(prob_samples, nextBest@overdose, bounds_closed = c(FALSE, TRUE)))
prob_target <- colMeans(h_in_range(prob_samples, nextBest@target))
# Eligible grid doses after accounting for maximum possible dose and discarding overdoses.
is_dose_eligible <- h_next_best_eligible_doses(data@doseGrid, doselimit, data@placebo, levels = TRUE) &
(prob_overdose <= nextBest@max_overdose_prob)
next_dose <- if (any(is_dose_eligible)) {
# If maximum target probability is higher than some numerical threshold,
# then take that level, otherwise stick to the maximum level that is OK.
# next_best_level is relative to eligible doses.
next_best_level <- ifelse(
test = any(prob_target[is_dose_eligible] > 0.05),
yes = which.max(prob_target[is_dose_eligible]),
no = sum(is_dose_eligible)
)
data@doseGrid[is_dose_eligible][next_best_level]
} else {
NA_real_
}
# Build plots, first for the target probability.
p1 <- ggplot() +
geom_bar(
data = data.frame(Dose = data@doseGrid, y = prob_target * 100),
aes(x = .data$Dose, y = .data$y),
stat = "identity",
position = "identity",
width = min(diff(data@doseGrid)) / 2,
colour = "darkgreen",
fill = "darkgreen"
) +
coord_cartesian(ylim = c(0, 100)) +
ylab(paste("Target probability [%]"))
if (is.finite(doselimit)) {
p1 <- p1 + geom_vline(xintercept = doselimit, lwd = 1.1, lty = 2, colour = "black")
}
if (any(is_dose_eligible)) {
p1 <- p1 +
geom_vline(xintercept = data@doseGrid[sum(is_dose_eligible)], lwd = 1.1, lty = 2, colour = "red") +
geom_point(
data = data.frame(x = next_dose, y = prob_target[is_dose_eligible][next_best_level] * 100 + 0.03),
aes(x = x, y = y),
size = 3,
pch = 25,
col = "red",
bg = "red"
)
}
# Second, for the overdosing probability.
p2 <- ggplot() +
geom_bar(
data = data.frame(Dose = data@doseGrid, y = prob_overdose * 100),
aes(x = .data$Dose, y = .data$y),
stat = "identity",
position = "identity",
width = min(diff(data@doseGrid)) / 2,
colour = "red",
fill = "red"
) +
geom_hline(yintercept = nextBest@max_overdose_prob * 100, lwd = 1.1, lty = 2, colour = "black") +
ylim(c(0, 100)) +
ylab("Overdose probability [%]")
# Place them below each other.
plot_joint <- gridExtra::arrangeGrob(p1, p2, nrow = 2)
list(
value = next_dose,
plot = plot_joint,
singlePlots = list(plot1 = p1, plot2 = p2),
probs = cbind(
dose = data@doseGrid,
target = prob_target,
overdose = prob_overdose
)
)
}
)
## NextBestNCRM-DataParts ----
#' @describeIn nextBest find the next best dose based on the NCRM method when
#' two parts trial is used.
#'
#' @aliases nextBest-NextBestNCRM-DataParts
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestNCRM-DataParts.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestNCRM",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "DataParts"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
# Exception when we are in part I or about to start part II!
if (all(data@part == 1L)) {
# Propose the highest possible dose (assuming that the dose limit came
# from reasonable increments rule, i.e. incrementsRelativeParts).
if (is.infinite(doselimit)) {
stop("A finite doselimit needs to be specified for Part I.")
}
list(value = doselimit, plot = NULL)
} else {
# Otherwise we will just do the standard thing.
callNextMethod(nextBest, doselimit, samples, model, data, ...)
}
}
)
## NextBestNCRMLoss ----
#' @describeIn nextBest find the next best dose based on the NCRM method and
#' loss function.
#'
#' @aliases nextBest-NextBestNCRMLoss
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestNCRMLoss.R
#'
setMethod("nextBest",
signature = signature(
nextBest = "NextBestNCRMLoss",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
# Matrix with samples from the dose-tox curve at the dose grid points.
prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
# Compute probabilities to be in target and overdose tox interval.
prob_underdosing <- colMeans(prob_samples < nextBest@target[1])
prob_target <- colMeans(h_in_range(prob_samples, nextBest@target))
prob_overdose <- colMeans(h_in_range(prob_samples, nextBest@overdose, bounds_closed = c(FALSE, TRUE)))
prob_mean <- colMeans(prob_samples)
prob_sd <- apply(prob_samples, 2, stats::sd)
is_unacceptable_specified <- any(nextBest@unacceptable != c(1, 1))
prob_mat <- if (!is_unacceptable_specified) {
cbind(underdosing = prob_underdosing, target = prob_target, overdose = prob_overdose)
} else {
prob_unacceptable <- colMeans(
h_in_range(prob_samples, nextBest@unacceptable, bounds_closed = c(FALSE, TRUE))
)
prob_excessive <- prob_overdose
prob_overdose <- prob_excessive + prob_unacceptable
cbind(
underdosing = prob_underdosing,
target = prob_target,
excessive = prob_excessive,
unacceptable = prob_unacceptable
)
}
posterior_loss <- as.vector(nextBest@losses %*% t(prob_mat))
names(posterior_loss) <- data@doseGrid
probs <- cbind(
dose = data@doseGrid,
prob_mat = prob_mat,
mean = prob_mean,
std_dev = prob_sd,
posterior_loss = posterior_loss
)
# Eligible grid doses after accounting for maximum possible dose and discarding overdoses.
is_dose_eligible <- h_next_best_eligible_doses(data@doseGrid, doselimit, data@placebo, levels = TRUE) &
(prob_overdose < nextBest@max_overdose_prob)
# Next best dose is the dose with the minimum loss function.
next_dose <- if (any(is_dose_eligible)) {
next_best_level <- which.min(posterior_loss[is_dose_eligible])
data@doseGrid[is_dose_eligible][next_best_level]
} else {
NA_real_
}
# Build plot.
p <- h_next_best_ncrm_loss_plot(
prob_mat = prob_mat,
posterior_loss = posterior_loss,
max_overdose_prob = nextBest@max_overdose_prob,
dose_grid = data@doseGrid,
max_eligible_dose_level = sum(is_dose_eligible),
doselimit = doselimit,
next_dose = next_dose,
is_unacceptable_specified = is_unacceptable_specified
)
c(list(value = next_dose, probs = probs), p)
}
)
## NextBestThreePlusThree ----
#' @describeIn nextBest find the next best dose based on the 3+3 method.
#'
#' @aliases nextBest-NextBestThreePlusThree
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestThreePlusThree.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestThreePlusThree",
doselimit = "missing",
samples = "missing",
model = "missing",
data = "Data"
),
definition = function(nextBest, doselimit, samples, model, data, ...) {
# The last dose level tested (not necessarily the maximum one).
last_level <- tail(data@xLevel, 1L)
# Get number of patients per grid's dose and DLT rate at the last level.
nPatients <- table(factor(data@x, levels = data@doseGrid))
n_dlts_last_level <- sum(data@y[data@xLevel == last_level])
dlt_rate_last_level <- n_dlts_last_level / nPatients[last_level]
level_change <- if (dlt_rate_last_level < 1 / 3) {
# Escalate it, unless this is the highest level or the higher dose was already tried.
ifelse((last_level == data@nGrid) || (nPatients[last_level + 1L] > 0), 0L, 1L)
} else {
# Rate is too high, deescalate it, unless an edge case of 1/3, where the decision
# depends on the num. of patients: if >3, then deescalate it, otherwise stay.
ifelse((dlt_rate_last_level == 1 / 3) && (nPatients[last_level] <= 3L), 0L, -1L)
}
next_dose_level <- last_level + level_change
# Do we stop here? Only if we have no MTD, or the next level has been tried
# enough (more than three patients already).
if (next_dose_level == 0L) {
next_dose <- NA
stop_here <- TRUE
} else {
next_dose <- data@doseGrid[next_dose_level]
stop_here <- nPatients[next_dose_level] > 3L
}
list(value = next_dose, stopHere = stop_here)
}
)
## NextBestDualEndpoint ----
#' @describeIn nextBest find the next best dose based on the dual endpoint
#' model. The additional list element `probs` contains the target and
#' overdosing probabilities (across all doses in the dose grid) used in the
#' derivation of the next best dose.
#'
#' @aliases nextBest-NextBestDualEndpoint
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestDualEndpoint.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestDualEndpoint",
doselimit = "numeric",
samples = "Samples",
model = "DualEndpoint",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
# Biomarker samples at the dose grid points.
biom_samples <- samples@data$betaW
prob_target <- if (nextBest@target_relative) {
# If 'Emax' parameter available, target biomarker level will be relative to 'Emax',
# otherwise, it will be relative to the maximum biomarker level achieved
# in dose range for a given sample.
if ("Emax" %in% names(samples)) {
lwr <- nextBest@target[1] * samples@data$Emax
upr <- nextBest@target[2] * samples@data$Emax
colMeans(apply(biom_samples, 2L, function(s) (s >= lwr) & (s <= upr)))
} else {
target_levels <- apply(biom_samples, 1L, function(x) {
rng <- range(x)
min(which(h_in_range(x, nextBest@target * diff(rng) + rng[1] + c(0, 1e-10), bounds_closed = c(FALSE, TRUE))))
})
prob_target <- as.vector(table(factor(target_levels, levels = 1:data@nGrid)))
prob_target / nrow(biom_samples)
}
} else {
colMeans(h_in_range(biom_samples, nextBest@target))
}
# Now, compute probabilities to be in overdose tox interval, then flag
# eligible grid doses after accounting for maximum possible dose and discarding overdoses.
prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples)
prob_overdose <- colMeans(h_in_range(prob_samples, nextBest@overdose, bounds_closed = c(FALSE, TRUE)))
is_dose_eligible <- h_next_best_eligible_doses(data@doseGrid, doselimit, data@placebo, levels = TRUE) &
(prob_overdose < nextBest@max_overdose_prob)
next_dose <- if (any(is_dose_eligible)) {
# If maximum target probability is higher the threshold, then take that
# level, otherwise stick to the maximum level that is eligible.
# next_dose_level is relative to eligible doses.
next_dose_level <- ifelse(
test = any(prob_target[is_dose_eligible] > nextBest@target_thresh),
yes = which.max(prob_target[is_dose_eligible]),
no = sum(is_dose_eligible)
)
data@doseGrid[is_dose_eligible][next_dose_level]
} else {
NA_real_
}
# Build plots, first for the target probability.
p1 <- ggplot() +
geom_bar(
data = data.frame(Dose = data@doseGrid, y = prob_target * 100),
aes(x = .data$Dose, y = .data$y),
stat = "identity",
position = "identity",
width = min(diff(data@doseGrid)) / 2,
colour = "darkgreen",
fill = "darkgreen"
) +
ylim(c(0, 100)) +
ylab(paste("Target probability [%]"))
if (is.finite(doselimit)) {
p1 <- p1 + geom_vline(xintercept = doselimit, lwd = 1.1, lty = 2, colour = "black")
}
if (any(is_dose_eligible)) {
p1 <- p1 +
geom_vline(xintercept = data@doseGrid[sum(is_dose_eligible)], lwd = 1.1, lty = 2, colour = "red") +
geom_point(
data = data.frame(x = next_dose, y = prob_target[is_dose_eligible][next_dose_level] * 100 + 0.03),
aes(x = x, y = y),
size = 3,
pch = 25,
col = "red",
bg = "red"
)
}
# Second, for the overdosing probability.
p2 <- ggplot() +
geom_bar(
data = data.frame(Dose = data@doseGrid, y = prob_overdose * 100),
aes(x = .data$Dose, y = .data$y),
stat = "identity",
position = "identity",
width = min(diff(data@doseGrid)) / 2,
colour = "red",
fill = "red"
) +
geom_hline(yintercept = nextBest@max_overdose_prob * 100, lwd = 1.1, lty = 2, colour = "black") +
ylim(c(0, 100)) +
ylab("Overdose probability [%]")
# Place them below each other.
plot_joint <- gridExtra::arrangeGrob(p1, p2, nrow = 2)
list(
value = next_dose,
plot = plot_joint,
singlePlots = list(plot1 = p1, plot2 = p2),
probs = cbind(dose = data@doseGrid, target = prob_target, overdose = prob_overdose)
)
}
)
## NextBestMinDist ----
#' @describeIn nextBest gives the dose which is below the dose limit and has an
#' estimated DLT probability which is closest to the target dose.
#'
#' @aliases nextBest-NextBestMinDist
#'
#' @export
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestMinDist",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
# Matrix with samples from the dose-tox curve at the dose grid points.
prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
dlt_prob <- colMeans(prob_samples)
# Determine the dose with the closest distance.
dose_target <- data@doseGrid[which.min(abs(dlt_prob - nextBest@target))]
# Determine next dose.
doses_eligible <- h_next_best_eligible_doses(
data@doseGrid,
doselimit,
data@placebo
)
next_dose_level_eligible <- which.min(abs(doses_eligible - dose_target))
next_dose <- doses_eligible[next_dose_level_eligible]
# Create a plot.
p <- ggplot(
data = data.frame(x = data@doseGrid, y = dlt_prob),
aes(.data$x, .data$y),
fill = "grey50",
colour = "grey50"
) +
geom_line() +
geom_point() +
coord_cartesian(xlim = range(data@doseGrid)) +
scale_x_continuous(
labels = as.character(data@doseGrid),
breaks = data@doseGrid,
guide = guide_axis(check.overlap = TRUE)
) +
geom_hline(yintercept = nextBest@target, linetype = "dotted") +
geom_vline(xintercept = dose_target, colour = "black", lwd = 1.1) +
geom_text(
data = data.frame(x = dose_target),
aes(.data$x, 0),
label = "Est",
vjust = -0.5,
hjust = 0.5,
colour = "black",
angle = 90
) +
xlab("Dose") +
ylab("Posterior toxicity probability")
if (is.finite(doselimit)) {
p <- p +
geom_vline(xintercept = doselimit, colour = "red", lwd = 1.1) +
geom_text(
data = data.frame(x = doselimit),
aes(.data$x, 0),
label = "Max",
vjust = -0.5,
hjust = -0.5,
colour = "red",
angle = 90
)
}
p <- p +
geom_vline(xintercept = next_dose, colour = "blue", lwd = 1.1) +
geom_text(
data = data.frame(x = next_dose),
aes(.data$x, 0),
label = "Next",
vjust = -0.5,
hjust = -1.5,
colour = "blue",
angle = 90
)
list(
value = next_dose,
probs = cbind(dose = data@doseGrid, dlt_prob = dlt_prob),
plot = p
)
}
)
## NextBestInfTheory ----
#' @describeIn nextBest gives the appropriate dose within an information
#' theoretic framework.
#'
#' @aliases nextBest-NextBestInfTheory
#'
#' @export
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestInfTheory",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
# Matrix with samples from the dose-tox curve at the dose grid points.
prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
criterion <- colMeans(h_info_theory_dist(prob_samples, nextBest@target, nextBest@asymmetry))
is_dose_eligible <- h_next_best_eligible_doses(data@doseGrid, doselimit, data@placebo, levels = TRUE)
doses_eligible <- data@doseGrid[is_dose_eligible]
next_best_level <- which.min(criterion[is_dose_eligible])
next_best <- doses_eligible[next_best_level]
list(value = next_best)
}
)
## NextBestTD ----
#' @describeIn nextBest find the next best dose based only on the DLT responses
#' and for [`LogisticIndepBeta`] model class object without DLT samples.
#'
#' @param model (`ModelTox`)\cr the DLT model.
#' @param in_sim (`flag`)\cr is this method used in simulations? Default as `FALSE`.
#' If this flag is `TRUE` and target dose estimates (during trial and end-of-trial)
#' are outside of the dose grid range, the information message is printed by
#' this method.
#'
#' @aliases nextBest-NextBestTD
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestTD.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestTD",
doselimit = "numeric",
samples = "missing",
model = "LogisticIndepBeta",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, model, data, in_sim = FALSE, ...) {
assert_flag(in_sim)
# 'drt' - during the trial, 'eot' end of trial.
prob_target_drt <- nextBest@prob_target_drt
prob_target_eot <- nextBest@prob_target_eot
# Target dose estimates, i.e. the dose with probability of the occurrence of
# a DLT that equals to the prob_target_drt or prob_target_eot.
dose_target_drt <- dose(x = prob_target_drt, model, ...)
dose_target_eot <- dose(x = prob_target_eot, model, ...)
# Find the next best doses in the doseGrid. The next best dose is the dose
# at level closest and below the target dose estimate.
# h_find_interval assumes that elements in doses_eligible are strictly increasing.
doses_eligible <- h_next_best_eligible_doses(data@doseGrid, doselimit, data@placebo)
next_dose_lev_drt <- h_find_interval(dose_target_drt, doses_eligible)
next_dose_drt <- doses_eligible[next_dose_lev_drt]
next_dose_lev_eot <- h_find_interval(dose_target_eot, doses_eligible)
next_dose_eot <- doses_eligible[next_dose_lev_eot]
# Find the variance of the log of the dose_target_eot.
mat <- matrix(
c(
-1 / (model@phi2),
-(log(prob_target_eot / (1 - prob_target_eot)) - model@phi1) / (model@phi2)^2
),
nrow = 1
)
var_dose_target_eot <- as.vector(mat %*% model@Pcov %*% t(mat))
# 95% credibility interval.
ci_dose_target_eot <- exp(log(dose_target_eot) + c(-1, 1) * 1.96 * sqrt(var_dose_target_eot))
cir_dose_target_eot <- ci_dose_target_eot[2] / ci_dose_target_eot[1]
# Build plot.
p <- h_next_best_td_plot(
prob_target_drt = prob_target_drt,
dose_target_drt = dose_target_drt,
prob_target_eot = prob_target_eot,
dose_target_eot = dose_target_eot,
data = data,
prob_dlt = prob(dose = data@doseGrid, model = model, ...),
doselimit = doselimit,
next_dose = next_dose_drt
)
if (!h_in_range(dose_target_drt, range = dose_grid_range(data), bounds_closed = TRUE) && !in_sim) {
warning(paste("TD", prob_target_drt * 100, "=", dose_target_drt, "not within dose grid"))
}
if (!h_in_range(dose_target_eot, range = dose_grid_range(data), bounds_closed = TRUE) && !in_sim) {
warning(paste("TD", prob_target_eot * 100, "=", dose_target_eot, "not within dose grid"))
}
list(
next_dose_drt = next_dose_drt,
prob_target_drt = prob_target_drt,
dose_target_drt = dose_target_drt,
next_dose_eot = next_dose_eot,
prob_target_eot = prob_target_eot,
dose_target_eot = dose_target_eot,
ci_dose_target_eot = ci_dose_target_eot,
ci_ratio_dose_target_eot = cir_dose_target_eot,
plot = p
)
}
)
## NextBestTDsamples ----
#' @describeIn nextBest find the next best dose based only on the DLT responses
#' and for [`LogisticIndepBeta`] model class object involving DLT samples.
#'
#' @aliases nextBest-NextBestTDsamples
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestTDsamples.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestTDsamples",
doselimit = "numeric",
samples = "Samples",
model = "LogisticIndepBeta",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, in_sim, ...) {
# Generate target dose samples, i.e. the doses with probability of the
# occurrence of a DLT that equals to the nextBest@prob_target_drt
# (or nextBest@prob_target_eot, respectively).
dose_target_drt_samples <- dose(x = nextBest@prob_target_drt, model, samples, ...)
dose_target_eot_samples <- dose(x = nextBest@prob_target_eot, model, samples, ...)
# Derive the prior/posterior estimates based on two above samples.
dose_target_drt <- nextBest@derive(dose_target_drt_samples)
dose_target_eot <- nextBest@derive(dose_target_eot_samples)
# Find the next doses in the doseGrid. The next dose is the dose at level
# closest and below the dose_target_drt (or dose_target_eot, respectively).
# h_find_interval assumes that elements in doses_eligible are strictly increasing.
doses_eligible <- h_next_best_eligible_doses(data@doseGrid, doselimit, data@placebo)
next_dose_lev_drt <- h_find_interval(dose_target_drt, doses_eligible)
next_dose_drt <- doses_eligible[next_dose_lev_drt]
next_dose_lev_eot <- h_find_interval(dose_target_eot, doses_eligible)
next_dose_eot <- doses_eligible[next_dose_lev_eot]
# 95% credibility interval.
ci_dose_target_eot <- as.numeric(quantile(dose_target_eot_samples, probs = c(0.025, 0.975)))
cir_dose_target_eot <- ci_dose_target_eot[2] / ci_dose_target_eot[1]
# Build plot.
p <- h_next_best_tdsamples_plot(
dose_target_drt_samples = dose_target_drt_samples,
dose_target_eot_samples = dose_target_eot_samples,
dose_target_drt = dose_target_drt,
dose_target_eot = dose_target_eot,
dose_grid_range = range(data@doseGrid),
nextBest = nextBest,
doselimit = doselimit,
next_dose = next_dose_drt
)
list(
next_dose_drt = next_dose_drt,
prob_target_drt = nextBest@prob_target_drt,
dose_target_drt = dose_target_drt,
next_dose_eot = next_dose_eot,
prob_target_eot = nextBest@prob_target_eot,
dose_target_eot = dose_target_eot,
ci_dose_target_eot = ci_dose_target_eot,
ci_ratio_dose_target_eot = cir_dose_target_eot,
plot = p
)
}
)
## NextBestMaxGain ----
#' @describeIn nextBest find the next best dose based only on pseudo DLT model
#' [`ModelTox`] and [`Effloglog`] efficacy model without samples.
#'
#' @param model (`ModelTox`)\cr the DLT model.
#' @param model_eff (`Effloglog`)\cr the efficacy model.
#' @param in_sim (`flag`)\cr is this method used in simulations? Default as `FALSE`.
#' If this flag is `TRUE` and target dose estimates (during trial and end-of-trial)
#' are outside of the dose grid range, the information message is printed by
#' this method.
#'
#' @aliases nextBest-NextBestMaxGain
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestMaxGain.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestMaxGain",
doselimit = "numeric",
samples = "missing",
model = "ModelTox",
data = "DataDual"
),
definition = function(nextBest, doselimit = Inf, model, data, model_eff, in_sim = FALSE, ...) {
assert_class(model_eff, "Effloglog")
assert_flag(in_sim)
# 'drt' - during the trial, 'eot' end of trial.
prob_target_drt <- nextBest@prob_target_drt
prob_target_eot <- nextBest@prob_target_eot
# Target dose estimates, i.e. the dose with probability of the occurrence of
# a DLT that equals to the prob_target_drt or prob_target_eot.
dose_target_drt <- dose(x = prob_target_drt, model, ...)
dose_target_eot <- dose(x = prob_target_eot, model, ...)
# Find the dose which gives the maximum gain.
dosegrid_range <- dose_grid_range(data)
opt <- optim(
par = dosegrid_range[1],
fn = function(DOSE) {
-gain(DOSE, model_dle = model, model_eff = model_eff, ...)
},
method = "L-BFGS-B",
lower = dosegrid_range[1],
upper = dosegrid_range[2]
)
dose_mg <- opt$par # this is G*. # no lintr
max_gain <- -opt$value
# Print info message if dose target is outside of the range.
if (!h_in_range(dose_target_drt, range = dose_grid_range(data), bounds_closed = FALSE) && !in_sim) {
print(paste("Estimated TD", prob_target_drt * 100, "=", dose_target_drt, "not within dose grid"))
}
if (!h_in_range(dose_target_eot, range = dose_grid_range(data), bounds_closed = FALSE) && !in_sim) {
print(paste("Estimated TD", prob_target_eot * 100, "=", dose_target_eot, "not within dose grid"))
}
if (!h_in_range(dose_mg, range = dose_grid_range(data), bounds_closed = FALSE) && !in_sim) {
print(paste("Estimated max gain dose =", dose_mg, "not within dose grid"))
}
# Get closest grid doses for a given target doses.
nb_doses_at_grid <- h_next_best_mg_doses_at_grid(
dose_target_drt = dose_target_drt,
dose_target_eot = dose_target_eot,
dose_mg = dose_mg,
dose_grid = data@doseGrid,
doselimit = doselimit,
placebo = data@placebo
)
# 95% credibility intervals and corresponding ratios for maximum gain dose and target dose eot.
ci <- h_next_best_mg_ci(
dose_target = dose_target_eot,
dose_mg = dose_mg,
prob_target = prob_target_eot,
placebo = data@placebo,
model = model,
model_eff = model_eff
)
# Build plot.
p <- h_next_best_mg_plot(
prob_target_drt = prob_target_drt,
dose_target_drt = dose_target_drt,
prob_target_eot = prob_target_eot,
dose_target_eot = dose_target_eot,
dose_mg = dose_mg,
max_gain = max_gain,
next_dose = nb_doses_at_grid$next_dose,
doselimit = doselimit,
data = data,
model = model,
model_eff = model_eff
)
list(
next_dose = nb_doses_at_grid$next_dose,
prob_target_drt = prob_target_drt,
dose_target_drt = dose_target_drt,
next_dose_drt = nb_doses_at_grid$next_dose_drt,
prob_target_eot = prob_target_eot,
dose_target_eot = dose_target_eot,
next_dose_eot = nb_doses_at_grid$next_dose_eot,
dose_max_gain = dose_mg,
next_dose_max_gain = nb_doses_at_grid$next_dose_mg,
ci_dose_target_eot = ci$ci_dose_target,
ci_ratio_dose_target_eot = ci$ci_ratio_dose_target,
ci_dose_max_gain = ci$ci_dose_mg,
ci_ratio_dose_max_gain = ci$ci_ratio_dose_mg,
plot = p
)
}
)
## NextBestMaxGainSamples ----
#' @describeIn nextBest find the next best dose based on DLT and efficacy
#' responses with DLT and efficacy samples.
#'
#' @param model (`ModelTox`)\cr the DLT model.
#' @param model_eff (`Effloglog` or `EffFlexi`)\cr the efficacy model.
#' @param samples_eff (`Samples`)\cr posterior samples from `model_eff` parameters
#' given `data`.
#' @param in_sim (`flag`)\cr is this method used in simulations? Default as `FALSE`.
#' If this flag is `TRUE` and target dose estimates (during trial and end-of-trial)
#' are outside of the dose grid range, the information message is printed by
#' this method.
#'
#' @aliases nextBest-NextBestMaxGainSamples
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestMaxGainSamples.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestMaxGainSamples",
doselimit = "numeric",
samples = "Samples",
model = "ModelTox",
data = "DataDual"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, model_eff, samples_eff, in_sim = FALSE, ...) {
assert_true(test_class(model_eff, "Effloglog") || test_class(model_eff, "EffFlexi"))
assert_class(samples_eff, "Samples")
assert_flag(in_sim)
# 'drt' - during the trial, 'eot' end of trial.
prob_target_drt <- nextBest@prob_target_drt
prob_target_eot <- nextBest@prob_target_eot
# Generate target dose samples, i.e. the doses with probability of the
# occurrence of a DLT that equals to the prob_target_drt or prob_target_eot.
dose_target_drt_samples <- dose(x = prob_target_drt, model, samples = samples, ...)
dose_target_eot_samples <- dose(x = prob_target_eot, model, samples = samples, ...)
# Derive the prior/posterior estimates based on two above samples.
dose_target_drt <- nextBest@derive(dose_target_drt_samples)
dose_target_eot <- nextBest@derive(dose_target_eot_samples)
# Gain samples.
gain_samples <- sapply(data@doseGrid, gain, model, samples, model_eff, samples_eff, ...)
# For every sample, get the dose (from the dose grid) that gives the maximum gain value.
dose_lev_mg_samples <- apply(gain_samples, 1, which.max)
dose_mg_samples <- data@doseGrid[dose_lev_mg_samples]
# Maximum gain dose estimate is the nth percentile of the maximum gain dose samples.
dose_mg <- nextBest@mg_derive(dose_mg_samples)
gain_values <- apply(gain_samples, 2, FUN = nextBest@mg_derive)
# Print info message if dose target is outside of the range.
dosegrid_range <- dose_grid_range(data)
if (!h_in_range(dose_target_drt, range = dosegrid_range, bounds_closed = FALSE) && !in_sim) {
print(paste("Estimated TD", prob_target_drt * 100, "=", dose_target_drt, "not within dose grid"))
}
if (!h_in_range(dose_target_eot, range = dosegrid_range, bounds_closed = FALSE) && !in_sim) {
print(paste("Estimated TD", prob_target_eot * 100, "=", dose_target_eot, "not within dose grid"))
}
if (!h_in_range(dose_mg, range = dosegrid_range, bounds_closed = FALSE) && !in_sim) {
print(paste("Estimated max gain dose =", dose_mg, "not within dose grid"))
}
# Get closest grid doses for a given target doses.
nb_doses_at_grid <- h_next_best_mg_doses_at_grid(
dose_target_drt = dose_target_drt,
dose_target_eot = dose_target_eot,
dose_mg = dose_mg,
dose_grid = data@doseGrid,
doselimit = doselimit,
placebo = data@placebo
)
# 95% credibility intervals and corresponding ratios for maximum gain dose and target dose eot.
ci_dose_mg <- as.numeric(quantile(dose_mg_samples, probs = c(0.025, 0.975)))
cir_dose_mg <- ci_dose_mg[2] / ci_dose_mg[1]
ci_dose_target_eot <- as.numeric(quantile(dose_target_eot, probs = c(0.025, 0.975)))
cir_dose_target_eot <- ci_dose_target_eot[2] / ci_dose_target_eot[1]
# Build plot.
p <- h_next_best_mgsamples_plot(
prob_target_drt = prob_target_drt,
dose_target_drt = dose_target_drt,
prob_target_eot = prob_target_eot,
dose_target_eot = dose_target_eot,
dose_mg = dose_mg,
dose_mg_samples = dose_mg_samples,
next_dose = nb_doses_at_grid$next_dose,
doselimit = doselimit,
dose_grid_range = dosegrid_range
)
list(
next_dose = nb_doses_at_grid$next_dose,
prob_target_drt = prob_target_drt,
dose_target_drt = dose_target_drt,
next_dose_drt = nb_doses_at_grid$next_dose_drt,
prob_target_eot = prob_target_eot,
dose_target_eot = dose_target_eot,
next_dose_eot = nb_doses_at_grid$next_dose_eot,
dose_max_gain = dose_mg,
next_dose_max_gain = nb_doses_at_grid$next_dose_mg,
ci_dose_target_eot = ci_dose_target_eot,
ci_ratio_dose_target_eot = cir_dose_target_eot,
ci_dose_max_gain = ci_dose_mg,
ci_ratio_dose_max_gain = cir_dose_mg,
plot = p
)
}
)
## NextBestProbMTDLTE ----
#' @describeIn nextBest find the next best dose based with the highest
#' probability of having a toxicity rate less or equal to the target toxicity
#' level.
#'
#' @aliases nextBest-NextBestProbMTDLTE
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestProbMTDLTE.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestProbMTDLTE",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit, samples, model, data, ...) {
# Matrix with samples from the dose-tox curve at the dose grid points.
prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
# Determine the maximum dose level with a toxicity probability below or
# equal to the target and calculate how often a dose is selected as MTD
# across iterations.
# The first element of the vector is the relative frequency that no
# dose in the grid is below or equal to the target, the
# second element that the 1st dose of the grid is the MTD, etc..
prob_mtd_lte <- prop.table(
table(factor(
rowSums(prob_samples <= nextBest@target),
levels = 0:data@nGrid
))
)
allocation_crit <- as.vector(prob_mtd_lte)
names(allocation_crit) <- as.character(c(0, data@doseGrid))
# In case that placebo is used, placebo and the portion that is not assigned
# to any dose of the grid are merged.
if (data@placebo) {
allocation_crit[1] <- sum(allocation_crit[1:2])
allocation_crit <- allocation_crit[-2]
}
# Handling of the portion that is not assigned to an active dose of
# the dose grid. The portion is added to the minimum active dose
# of the dose grid.
allocation_crit[2] <- sum(allocation_crit[1:2])
allocation_crit <- allocation_crit[-1]
# Determine the dose with the highest relative frequency.
allocation_crit_dose <- as.numeric(names(allocation_crit))
dose_target <- allocation_crit_dose[which.max(allocation_crit)]
# Determine next dose.
doses_eligible <- h_next_best_eligible_doses(
data@doseGrid,
doselimit,
data@placebo
)
next_dose_level_eligible <- which.min(abs(doses_eligible - dose_target))
next_dose <- doses_eligible[next_dose_level_eligible]
# Create a plot.
plt_data <- if (data@placebo && (data@doseGrid[1] == next_dose)) {
data.frame(
x = as.factor(data@doseGrid),
y = c(0, as.numeric(allocation_crit)) * 100
)
} else {
data.frame(
x = as.factor(allocation_crit_dose),
y = as.numeric(allocation_crit) * 100
)
}
p <- ggplot(
data = plt_data,
fill = "grey50",
colour = "grey50"
) +
geom_col(aes(x, y), fill = "grey75") +
scale_x_discrete(drop = FALSE, guide = guide_axis(check.overlap = TRUE)) +
geom_vline(
xintercept = as.factor(dose_target),
lwd = 1.1,
colour = "black"
) +
geom_text(
data = data.frame(x = as.factor(dose_target)),
aes(.data$x, 0),
label = "Est",
vjust = -0.5,
hjust = -0.5,
colour = "black",
angle = 90
) +
xlab("Dose") +
ylab(paste("Allocation criterion [%]"))
if (is.finite(doselimit)) {
doselimit_level <- if (sum(allocation_crit_dose == doselimit) > 0) {
which(allocation_crit_dose == doselimit)
} else {
ifelse(test = data@placebo && (data@doseGrid[1] == next_dose),
yes = 1.5,
no = sum(allocation_crit_dose < doselimit) + 0.5
)
}
p <- p +
geom_vline(
xintercept = doselimit_level,
colour = "red", lwd = 1.1
) +
geom_text(
data = data.frame(x = doselimit_level),
aes(.data$x, 0),
label = "Max",
vjust = -0.5,
hjust = -1.5,
colour = "red",
angle = 90
)
}
p <- p +
geom_vline(
xintercept = as.factor(next_dose),
colour = "blue", lwd = 1.1
) +
geom_text(
data = data.frame(x = as.factor(next_dose)),
aes(.data$x, 0),
label = "Next",
vjust = -0.5,
hjust = -2.5,
colour = "blue",
angle = 90
)
list(
value = next_dose,
allocation = cbind(dose = allocation_crit_dose, allocation = allocation_crit),
plot = p
)
}
)
## NextBestProbMTDMinDist ----
#' @describeIn nextBest find the next best dose based with the highest
#' probability of having a toxicity rate with minimum distance to the
#' target toxicity level.
#'
#' @aliases nextBest-NextBestProbMTDMinDist
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestProbMtdMinDist.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestProbMTDMinDist",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit, samples, model, data, ...) {
# Matrix with samples from the dose-tox curve at the dose grid points.
prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
# Determine which dose level has the minimum distance to target.
dose_min_mtd_dist <- apply(
prob_samples, 1, function(x) which.min(abs(x - nextBest@target))
)
allocation_crit <- prop.table(
table(factor(dose_min_mtd_dist, levels = 1:data@nGrid))
)
names(allocation_crit) <- as.character(data@doseGrid)
# In case that placebo is used, placebo and the first non-placebo dose
# of the grid are merged.
if (data@placebo) {
allocation_crit[2] <- sum(allocation_crit[1:2])
allocation_crit <- allocation_crit[-1]
}
# Determine the dose with the highest relative frequency.
allocation_crit_dose <- as.numeric(names(allocation_crit))
dose_target <- allocation_crit_dose[which.max(allocation_crit)]
# Determine next dose.
doses_eligible <- h_next_best_eligible_doses(
data@doseGrid,
doselimit,
data@placebo
)
next_dose_level_eligible <- which.min(abs(doses_eligible - dose_target))
next_dose <- doses_eligible[next_dose_level_eligible]
# Create a plot.
plt_data <- if (data@placebo && data@doseGrid[1] == next_dose) {
data.frame(
x = as.factor(data@doseGrid),
y = c(0, as.numeric(allocation_crit)) * 100
)
} else {
data.frame(
x = as.factor(allocation_crit_dose),
y = as.numeric(allocation_crit) * 100
)
}
p <- ggplot(
data = plt_data,
fill = "grey50",
colour = "grey50"
) +
geom_col(aes(x, y), fill = "grey75") +
scale_x_discrete(guide = guide_axis(check.overlap = TRUE)) +
geom_vline(
xintercept = as.factor(dose_target),
lwd = 1.1,
colour = "black"
) +
geom_text(
data = data.frame(x = as.factor(dose_target)),
aes(.data$x, 0),
label = "Est",
vjust = -0.5,
hjust = -0.5,
colour = "black",
angle = 90
) +
xlab("Dose") +
ylab(paste("Allocation criterion [%]"))
if (is.finite(doselimit)) {
doselimit_level <- if (any(allocation_crit_dose == doselimit)) {
which(allocation_crit_dose == doselimit)
} else {
ifelse(test = data@placebo && data@doseGrid[1] == next_dose,
yes = 1.5,
no = sum(allocation_crit_dose < doselimit) + 0.5
)
}
p <- p +
geom_vline(
xintercept = doselimit_level,
colour = "red", lwd = 1.1
) +
geom_text(
data = data.frame(x = doselimit_level),
aes(.data$x, 0),
label = "Max",
vjust = -0.5,
hjust = -1.5,
colour = "red",
angle = 90
)
}
p <- p +
geom_vline(
xintercept = as.factor(next_dose),
colour = "blue", lwd = 1.1
) +
geom_text(
data = data.frame(x = as.factor(next_dose)),
aes(.data$x, 0),
label = "Next",
vjust = -0.5,
hjust = -2.5,
colour = "blue",
angle = 90
)
list(
value = next_dose,
allocation = cbind(dose = allocation_crit_dose, allocation = allocation_crit),
plot = p
)
}
)
## NextBestOrdinal ----
#' @describeIn nextBest find the next best dose for ordinal CRM models.
#'
#' @aliases nextBest-NextBestOrdinal
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestOrdinal.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestOrdinal",
doselimit = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "Data"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
stop(
paste0(
"NextBestOrdinal objects can only be used with LogisticLogNormalOrdinal ",
"models and DataOrdinal data objects. In this case, the model is a '",
class(model),
"' object and the data is in a ",
class(data),
" object."
)
)
}
)
#' @describeIn nextBest find the next best dose for ordinal CRM models.
#'
#' @aliases nextBest-NextBestOrdinal
#'
#' @export
#' @example examples/Rules-method-nextBest-NextBestOrdinal.R
#'
setMethod(
f = "nextBest",
signature = signature(
nextBest = "NextBestOrdinal",
doselimit = "numeric",
samples = "Samples",
model = "LogisticLogNormalOrdinal",
data = "DataOrdinal"
),
definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
nextBest(
nextBest = nextBest@rule,
doselimit = doselimit,
samples = h_convert_ordinal_samples(samples, nextBest@grade),
model = h_convert_ordinal_model(model, nextBest@grade),
data = h_convert_ordinal_data(data, nextBest@grade),
...
)
}
)
# maxDose ----
## generic ----
#' Determine the Maximum Possible Next Dose
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This function determines the upper limit of the next dose based on the
#' `increments`and the `data`.
#'
#' @param increments (`Increments`)\cr the rule for the next best dose.
#' @param data (`Data`)\cr input data.
#' @param ... additional arguments without method dispatch.
#'
#' @return A `number`, the maximum possible next dose.
#'
#' @export
#'
setGeneric(
name = "maxDose",
def = function(increments, data, ...) {
standardGeneric("maxDose")
},
valueClass = "numeric"
)
## IncrementsRelative ----
#' @describeIn maxDose determine the maximum possible next dose based on
#' relative increments.
#'
#' @aliases maxDose-IncrementsRelative
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsRelative.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsRelative",
data = "Data"
),
definition = function(increments, data, ...) {
last_dose <- data@x[data@nObs]
# Determine in which interval the `last_dose` is.
assert_true(last_dose >= head(increments@intervals, 1))
last_dose_interval <- findInterval(x = last_dose, vec = increments@intervals)
(1 + increments@increments[last_dose_interval]) * last_dose
}
)
## IncrementsRelativeDLT ----
#' @describeIn maxDose determine the maximum possible next dose based on
#' relative increments determined by DLTs so far.
#'
#' @aliases maxDose-IncrementsRelativeDLT
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsRelativeDLT.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsRelativeDLT",
data = "Data"
),
definition = function(increments, data, ...) {
dlt_count <- sum(data@y)
# Determine in which interval the `dlt_count` is.
assert_true(dlt_count >= increments@intervals[1])
dlt_count_interval <- findInterval(x = dlt_count, vec = increments@intervals)
(1 + increments@increments[dlt_count_interval]) * data@x[data@nObs]
}
)
## IncrementsRelativeDLTCurrent ----
#' @describeIn maxDose determine the maximum possible next dose based on
#' relative increments determined by DLTs in the current cohort.
#'
#' @aliases maxDose-IncrementsRelativeDLTCurrent
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsRelativeDLTCurrent.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsRelativeDLTCurrent",
data = "Data"
),
definition = function(increments, data, ...) {
last_dose <- data@x[data@nObs]
# Determine how many DLTs have occurred in the last cohort.
last_cohort <- data@cohort[data@nObs]
last_cohort_indices <- which(data@cohort == last_cohort)
dlt_count_lcohort <- sum(data@y[last_cohort_indices])
# Determine in which interval the `dlt_count_lcohort` is.
assert_true(dlt_count_lcohort >= increments@intervals[1])
dlt_count_lcohort_int <- findInterval(x = dlt_count_lcohort, vec = increments@intervals)
(1 + increments@increments[dlt_count_lcohort_int]) * last_dose
}
)
## IncrementsRelativeParts ----
#' @describeIn maxDose determine the maximum possible next dose based on
#' relative increments as well as part 1 and beginning of part 2.
#'
#' @aliases maxDose-IncrementsRelativeParts
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsRelativeParts.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsRelativeParts",
data = "DataParts"
),
definition = function(increments, data, ...) {
all_in_part1 <- all(data@part == 1L)
incrmnt <- if (all_in_part1) {
part2_started <- data@nextPart == 2L
if (part2_started) {
any_dlt <- any(data@y == 1L)
if (any_dlt) {
increments@dlt_start
} else if (increments@clean_start <= 0L) {
increments@clean_start
}
} else {
1L
}
}
if (is.null(incrmnt)) {
callNextMethod(increments, data, ...)
} else {
max_dose_lev_part1 <- match_within_tolerance(max(data@x), data@part1Ladder)
new_max_dose_level <- max_dose_lev_part1 + incrmnt
assert_true(new_max_dose_level >= 0L)
assert_true(new_max_dose_level <= length(data@part1Ladder))
data@part1Ladder[new_max_dose_level]
}
}
)
## IncrementsDoseLevels ----
#' @describeIn maxDose determine the maximum possible next dose based on
#' the number of dose grid levels. That is, the max dose is determined as
#' the one which level is equal to: base dose level + level increment.
#' The base dose level is the level of the last dose in grid or the level
#' of the maximum dose applied, which is defined in `increments` object.
#' Find out more in [`IncrementsDoseLevels`].
#'
#' @aliases maxDose-IncrementsDoseLevels
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsDoseLevels.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsDoseLevels",
data = "Data"
),
definition = function(increments, data, ...) {
# Determine what is the basis level for increment,
# i.e. the last dose or the max dose applied.
basis_dose_level <- ifelse(
increments@basis_level == "last", data@xLevel[data@nObs], max(data@xLevel)
)
max_dose_level <- min(basis_dose_level + increments@levels, data@nGrid)
data@doseGrid[max_dose_level]
}
)
## IncrementsHSRBeta ----
#' @describeIn maxDose determine the maximum possible next dose for escalation.
#'
#' @aliases maxDose-IncrementsHSRBeta
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsHSRBeta.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsHSRBeta",
data = "Data"
),
definition = function(increments, data, ...) {
# Summary of observed data per dose level.
y <- factor(data@y, levels = c("0", "1"))
dlt_tab <- table(y, data@x)
# Ignore placebo if applied.
if (data@placebo == TRUE & min(data@x) == data@doseGrid[1]) {
dlt_tab <- dlt_tab[, -1]
}
# Extract dose names as these get lost if only one dose available.
non_plcb_doses <- unique(sort(as.numeric(colnames(dlt_tab))))
# Toxicity probability per dose level.
x <- dlt_tab[2, ]
n <- apply(dlt_tab, 2, sum)
tox_prob <- pbeta(
increments@target,
x + increments@a,
n - x + increments@b,
lower.tail = FALSE
)
# Return the min toxic dose level or maximum dose level if no dose is toxic,
# while ignoring placebo.
dose_tox <- if (sum(tox_prob > increments@prob) > 0) {
min(non_plcb_doses[which(tox_prob > increments@prob)])
} else {
# Add small value to max dose, so that the max dose is always smaller.
max(data@doseGrid) + 0.01
}
# Determine the next maximum possible dose.
# In case that the first active dose is above probability threshold,
# the first active dose is reported as maximum. I.e. in case that placebo is used,
# the second dose is reported. Please note that this rule should be used together
# with the hard safety stopping rule to avoid inconsistent results.
max(data@doseGrid[data@doseGrid < dose_tox], data@doseGrid[data@placebo + 1])
}
)
## IncrementsMin ----
#' @describeIn maxDose determine the maximum possible next dose based on
#' multiple increment rules, taking the minimum across individual increments.
#'
#' @aliases maxDose-IncrementsMin
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsMin.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsMin",
data = "Data"
),
definition = function(increments, data, ...) {
individual_results <- sapply(increments@increments_list, maxDose, data = data, ...)
min(individual_results)
}
)
#' @describeIn maxDose determine the maximum possible next dose based on
#' multiple increment rules, taking the minimum across individual increments.
#'
#' @aliases maxDose-IncrementsMin
#'
#' @export
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsMin",
data = "DataOrdinal"
),
definition = function(increments, data, ...) {
individual_results <- sapply(increments@increments_list, maxDose, data = data, ...)
min(individual_results)
}
)
## IncrementsOrdinal ----
#' @describeIn maxDose determine the maximum possible next dose in an ordinal
#' CRM trial
#'
#' @aliases maxDose-IncrementsOrdinal
#'
#' @export
#' @example examples/Rules-method-maxDose-IncrementsOrdinal.R
#'
setMethod(
f = "maxDose",
signature = signature(
increments = "IncrementsOrdinal",
data = "DataOrdinal"
),
definition = function(increments, data, ...) {
maxDose(
increments = increments@rule,
data = h_convert_ordinal_data(
data,
increments@grade,
...
)
)
}
)
# nolint start
## ============================================================
## --------------------------------------------------
## "AND" combination of stopping rules
## --------------------------------------------------
##' The method combining two atomic stopping rules
##'
##' @param e1 First \code{\linkS4class{Stopping}} object
##' @param e2 Second \code{\linkS4class{Stopping}} object
##' @return The \code{\linkS4class{StoppingAll}} object
##'
##' @example examples/Rules-method-and-stopping-stopping.R
##' @keywords methods
setMethod("&",
signature(
e1 = "Stopping",
e2 = "Stopping"
),
def =
function(e1, e2) {
StoppingAll(list(e1, e2))
}
)
##' The method combining a stopping list and an atomic
##'
##' @param e1 \code{\linkS4class{StoppingAll}} object
##' @param e2 \code{\linkS4class{Stopping}} object
##' @return The modified \code{\linkS4class{StoppingAll}} object
##'
##' @example examples/Rules-method-and-stoppingAll-stopping.R
##' @keywords methods
setMethod("&",
signature(
e1 = "StoppingAll",
e2 = "Stopping"
),
def =
function(e1, e2) {
e1@stop_list <- c(
e1@stop_list,
e2
)
return(e1)
}
)
##' The method combining an atomic and a stopping list
##'
##' @param e1 \code{\linkS4class{Stopping}} object
##' @param e2 \code{\linkS4class{StoppingAll}} object
##' @return The modified \code{\linkS4class{StoppingAll}} object
##'
##' @example examples/Rules-method-and-stopping-stoppingAll.R
##' @keywords methods
setMethod("&",
signature(
e1 = "Stopping",
e2 = "StoppingAll"
),
def =
function(e1, e2) {
e2@stop_list <- c(
e1,
e2@stop_list
)
return(e2)
}
)
## --------------------------------------------------
## "OR" combination of stopping rules
## --------------------------------------------------
##' The method combining two atomic stopping rules
##'
##' @param e1 First \code{\linkS4class{Stopping}} object
##' @param e2 Second \code{\linkS4class{Stopping}} object
##' @return The \code{\linkS4class{StoppingAny}} object
##'
##' @aliases |,Stopping,Stopping-method
##' @name or-Stopping-Stopping
##' @example examples/Rules-method-or-stopping-stopping.R
##' @keywords methods
setMethod("|",
signature(
e1 = "Stopping",
e2 = "Stopping"
),
def =
function(e1, e2) {
StoppingAny(list(e1, e2))
}
)
##' The method combining a stopping list and an atomic
##'
##' @param e1 \code{\linkS4class{StoppingAny}} object
##' @param e2 \code{\linkS4class{Stopping}} object
##' @return The modified \code{\linkS4class{StoppingAny}} object
##'
##' @aliases |,StoppingAny,Stopping-method
##' @name or-Stopping-StoppingAny
##' @example examples/Rules-method-or-stoppingAny-stopping.R
##' @keywords methods
setMethod("|",
signature(
e1 = "StoppingAny",
e2 = "Stopping"
),
def =
function(e1, e2) {
e1@stop_list <- c(
e1@stop_list,
e2
)
return(e1)
}
)
##' The method combining an atomic and a stopping list
##'
##' @param e1 \code{\linkS4class{Stopping}} object
##' @param e2 \code{\linkS4class{StoppingAny}} object
##' @return The modified \code{\linkS4class{StoppingAny}} object
##'
##' @aliases |,Stopping,StoppingAny-method
##' @name or-StoppingAny-Stopping
##' @example examples/Rules-method-or-stopping-stoppingAny.R
##' @keywords methods
setMethod("|",
signature(
e1 = "Stopping",
e2 = "StoppingAny"
),
def =
function(e1, e2) {
e2@stop_list <- c(
e1,
e2@stop_list
)
return(e2)
}
)
# nolint end
# Stopping ----
## generic ----
#' Stop the trial?
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This function returns whether to stop the trial.
#'
#' @param stopping (`Stopping`)\cr the rule for stopping the trial.
#' @param dose the recommended next best dose.
#' @param samples (`Samples`)\cr the mcmc samples.
#' @param model (`GeneralModel`)\cr the model.
#' @param data (`Data`)\cr input data.
#' @param ... additional arguments without method dispatch.
#'
#' @return logical value: `TRUE` if the trial can be stopped, `FALSE`
#' otherwise. It should have an attribute `message` which gives the reason
#' for the decision.
#'
#' @export
#' @example examples/Rules-method-CombiningStoppingRulesAndOr.R
setGeneric(
name = "stopTrial",
def = function(stopping, dose, samples, model, data, ...) {
standardGeneric("stopTrial")
},
valueClass = "logical"
)
## StoppingMissingDose ----
#' @describeIn stopTrial Stop based on value returned by next best dose.
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' @aliases stopTrial-StoppingMissingDose
#' @example examples/Rules-method-stopTrial-StoppingMissingDose.R
#'
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingMissingDose",
dose = "numeric",
samples = "ANY",
model = "ANY",
data = "Data"
),
definition = function(stopping, dose, samples, model, data, ...) {
do_stop <- is.na(dose) || (data@placebo && dose == min(data@doseGrid))
msg <- paste(
"Next dose is",
ifelse(
do_stop,
paste(
ifelse(
data@placebo && dose == min(data@doseGrid),
"placebo dose",
"NA"
),
", i.e., no active dose is safe enough according to the NextBest rule."
),
"available at the dose grid."
)
)
structure(do_stop,
message = msg,
report_label = stopping@report_label
)
}
)
# nolint start
## --------------------------------------------------
## Stopping based on multiple stopping rules
## --------------------------------------------------
##' @describeIn stopTrial Stop based on multiple stopping rules
##' @example examples/Rules-method-stopTrial-StoppingList.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingList",
dose = "ANY",
samples = "ANY",
model = "ANY",
data = "ANY"
),
def =
function(stopping, dose, samples, model, data, ...) {
## evaluate the individual stopping rules
## in the list
individualResults <-
if (missing(samples)) {
lapply(stopping@stop_list,
stopTrial,
dose = dose,
model = model,
data = data,
...
)
} else {
lapply(stopping@stop_list,
stopTrial,
dose = dose,
samples = samples,
model = model,
data = data,
...
)
}
## summarize to obtain overall result
overallResult <- stopping@summary(as.logical(individualResults))
## retrieve individual text messages,
## but let them in the list structure
overallText <- lapply(individualResults, attr, "message")
return(structure(overallResult,
message = overallText,
individual = individualResults
))
}
)
## --------------------------------------------------
## Stopping based on fulfillment of all multiple stopping rules
## --------------------------------------------------
##' @describeIn stopTrial Stop based on fulfillment of all multiple stopping
##' rules
##'
##' @example examples/Rules-method-stopTrial-StoppingAll.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingAll",
dose = "ANY",
samples = "ANY",
model = "ANY",
data = "ANY"
),
def =
function(stopping, dose, samples, model, data, ...) {
## evaluate the individual stopping rules
## in the list
individualResults <-
if (missing(samples)) {
lapply(stopping@stop_list,
stopTrial,
dose = dose,
model = model,
data = data,
...
)
} else {
lapply(stopping@stop_list,
stopTrial,
dose = dose,
samples = samples,
model = model,
data = data,
...
)
}
## summarize to obtain overall result
overallResult <- all(as.logical(individualResults))
## retrieve individual text messages,
## but let them in the list structure
overallText <- lapply(individualResults, attr, "message")
return(structure(overallResult,
message = overallText,
individual = individualResults,
report_label = stopping@report_label
))
}
)
## --------------------------------------------------
## Stopping based on fulfillment of any stopping rule
## --------------------------------------------------
##' @describeIn stopTrial Stop based on fulfillment of any stopping rule
##'
##' @example examples/Rules-method-stopTrial-StoppingAny.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingAny",
dose = "ANY",
samples = "ANY",
model = "ANY",
data = "ANY"
),
def =
function(stopping, dose, samples, model, data, ...) {
## evaluate the individual stopping rules
## in the list
individualResults <-
if (missing(samples)) {
lapply(stopping@stop_list,
stopTrial,
dose = dose,
model = model,
data = data,
...
)
} else {
lapply(stopping@stop_list,
stopTrial,
dose = dose,
samples = samples,
model = model,
data = data,
...
)
}
## summarize to obtain overall result
overallResult <- any(as.logical(individualResults))
## retrieve individual text messages,
## but let them in the list structure
overallText <- lapply(individualResults, attr, "message")
return(structure(overallResult,
message = overallText,
individual = individualResults,
report_label = stopping@report_label
))
}
)
## --------------------------------------------------
## Stopping based on number of cohorts near to next best dose
## --------------------------------------------------
##' @describeIn stopTrial Stop based on number of cohorts near to next best dose
##'
##' @example examples/Rules-method-stopTrial-StoppingCohortsNearDose.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingCohortsNearDose",
dose = "numeric",
samples = "ANY",
model = "ANY",
data = "Data"
),
def =
function(stopping, dose, samples, model, data, ...) {
## determine the range where the cohorts must lie in
lower <- (100 - stopping@percentage) / 100 * dose
upper <- (100 + stopping@percentage) / 100 * dose
## which patients lie there?
indexPatients <- which((data@x >= lower) & (data@x <= upper))
## how many cohorts?
nCohorts <- length(unique(data@cohort[indexPatients]))
## so can we stop?
doStop <- nCohorts >= stopping@nCohorts
## generate message
text <- paste(nCohorts,
" cohorts lie within ",
stopping@percentage,
"% of the next best dose ",
dose,
". This ",
ifelse(doStop, "reached", "is below"),
" the required ",
stopping@nCohorts,
" cohorts",
sep = ""
)
## return both
return(structure(doStop,
message = text,
report_label = stopping@report_label
))
}
)
## -------------------------------------------------------------
## Stopping based on number of patients near to next best dose
## -------------------------------------------------------------
##' @describeIn stopTrial Stop based on number of patients near to next best
##' dose
##'
##' @example examples/Rules-method-stopTrial-StoppingPatientsNearDose.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingPatientsNearDose",
dose = "numeric",
samples = "ANY",
model = "ANY",
data = "Data"
),
def =
function(stopping, dose, samples, model, data, ...) {
## determine the range where the cohorts must lie in
lower <- (100 - stopping@percentage) / 100 * dose
upper <- (100 + stopping@percentage) / 100 * dose
## how many patients lie there?
nPatients <- ifelse(
is.na(dose),
0,
sum((data@x >= lower) & (data@x <= upper))
)
## so can we stop?
doStop <- nPatients >= stopping@nPatients
## generate message
text <- paste(nPatients,
" patients lie within ",
stopping@percentage,
"% of the next best dose ",
dose,
". This ",
ifelse(doStop, "reached", "is below"),
" the required ",
stopping@nPatients,
" patients",
sep = ""
)
## return both
return(structure(doStop,
message = text,
report_label = stopping@report_label
))
}
)
## --------------------------------------------------
## Stopping based on minimum number of cohorts
## --------------------------------------------------
##' @describeIn stopTrial Stop based on minimum number of cohorts
##'
##' @example examples/Rules-method-stopTrial-StoppingMinCohorts.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingMinCohorts",
dose = "ANY",
samples = "ANY",
model = "ANY",
data = "Data"
),
def =
function(stopping, dose, samples, model, data, ...) {
## determine number of cohorts
nCohorts <- length(unique(data@cohort))
## so can we stop?
doStop <- nCohorts >= stopping@nCohorts
## generate message
text <-
paste(
"Number of cohorts is",
nCohorts,
"and thus",
ifelse(doStop, "reached", "below"),
"the prespecified minimum number",
stopping@nCohorts
)
## return both
return(structure(doStop,
message = text,
report_label = stopping@report_label
))
}
)
## --------------------------------------------------
## Stopping based on minimum number of patients
## --------------------------------------------------
##' @describeIn stopTrial Stop based on minimum number of patients
##'
##' @example examples/Rules-method-stopTrial-StoppingMinPatients.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingMinPatients",
dose = "ANY",
samples = "ANY",
model = "ANY",
data = "Data"
),
def =
function(stopping, dose, samples, model, data, ...) {
## so can we stop?
doStop <- data@nObs >= stopping@nPatients
## generate message
text <-
paste(
"Number of patients is",
data@nObs,
"and thus",
ifelse(doStop, "reached", "below"),
"the prespecified minimum number",
stopping@nPatients
)
## return both
return(structure(doStop,
message = text,
report_label = stopping@report_label
))
}
)
# nolint end
## StoppingTargetProb ----
#' @describeIn stopTrial Stop based on probability of target tox interval
#'
#' @aliases stopTrial-StoppingTargetProb
#' @example examples/Rules-method-stopTrial-StoppingTargetProb.R
setMethod(
f = "stopTrial",
signature =
signature(
stopping = "StoppingTargetProb",
dose = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "ANY"
),
definition = function(stopping, dose, samples, model, data, ...) {
# Compute probability to be in target interval.
prob_target <- ifelse(
is.na(dose),
0,
mean(
prob(dose = dose, model, samples, ...) >= stopping@target[1] &
prob(dose = dose, model, samples, ...) <= stopping@target[2]
)
)
do_stop <- prob_target >= stopping@prob
msg <- paste(
"Probability for target toxicity is",
round(prob_target * 100),
"% for dose",
dose,
"and thus",
ifelse(do_stop, "above", "below"),
"the required",
round(stopping@prob * 100),
"%"
)
structure(
do_stop,
message = msg,
report_label = stopping@report_label
)
}
)
# nolint start
## --------------------------------------------------
## Stopping based on MTD distribution
## --------------------------------------------------
##' @describeIn stopTrial Stop based on MTD distribution
##'
##' @example examples/Rules-method-stopTrial-StoppingMTDdistribution.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingMTDdistribution",
dose = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "ANY"
),
def =
function(stopping, dose, samples, model, data, ...) {
## First, generate the MTD samples.
## add prior data and samples to the
## function environment so that they
## can be used.
mtdSamples <- dose(
x = stopping@target,
model,
samples,
...
)
## what is the absolute threshold?
absThresh <- stopping@thresh * dose
## what is the probability to be above this dose?
prob <- ifelse(
is.na(absThresh),
0,
mean(mtdSamples > absThresh)
)
## so can we stop?
doStop <- prob >= stopping@prob
## generate message
text <-
paste(
"Probability of MTD above",
round(stopping@thresh * 100),
"% of current dose",
dose,
"is",
round(prob * 100),
"% and thus",
ifelse(doStop, "greater than or equal to", "strictly less than"),
"the required",
round(stopping@prob * 100),
"%"
)
## return both
return(structure(doStop,
message = text,
report_label = stopping@report_label
))
}
)
# nolint end
## StoppingMTDCV ----
#' @rdname stopTrial
#'
#' @description Stopping rule based precision of the MTD estimation.
#' The trial is stopped, when the MTD can be estimated with sufficient precision.
#' The criteria is based on the robust coefficient of variation (CV) calculated
#' from the posterior distribution.
#' The robust CV is defined `mad(MTD) / median(MTD)`, where `mad` is the median
#' absolute deviation.
#'
#' @aliases stopTrial-StoppingMTDCV
#' @example examples/Rules-method-stopTrial-StoppingMTDCV.R
#' @export
#'
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingMTDCV",
dose = "numeric",
samples = "Samples",
model = "GeneralModel",
data = "ANY"
),
definition = function(stopping, dose, samples, model, data, ...) {
mtd_samples <- dose(
x = stopping@target,
model,
samples,
...
)
# CV of MTD expressed as percentage, derived based on MTD posterior samples.
mtd_cv <- (mad(mtd_samples) / median(mtd_samples)) * 100
do_stop <- mtd_cv <= stopping@thresh_cv
msg <- paste(
"CV of MTD is",
round(mtd_cv),
"% and thus",
ifelse(do_stop, "below", "above"),
"the required precision threshold of",
round(stopping@thresh_cv),
"%"
)
structure(
do_stop,
message = msg,
report_label = stopping@report_label
)
}
)
## StoppingLowestDoseHSRBeta ----
#' @rdname stopTrial
#'
#' @description Stopping based based on the lowest non placebo dose. The trial is
#' stopped when the lowest non placebo dose meets the Hard
#' Safety Rule, i.e. it is deemed to be overly toxic. Stopping is based on the
#' observed data at the lowest dose level using a Bin-Beta model
#' based on DLT probability.
#'
#' @aliases stopTrial-StoppingLowestDoseHSRBeta
#' @example examples/Rules-method-stopTrial-StoppingLowestDoseHSRBeta.R
#' @export
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingLowestDoseHSRBeta",
dose = "numeric",
samples = "Samples"
),
definition = function(stopping, dose, samples, model, data, ...) {
# Actual number of patients at first active dose.
n <- sum(data@x == data@doseGrid[data@placebo + 1])
# Determine toxicity probability of the first active dose.
tox_prob_first_dose <-
if (n > 0) {
x <- sum(data@y[which(data@x == data@doseGrid[data@placebo + 1])])
pbeta(stopping@target, x + stopping@a, n - x + stopping@b, lower.tail = FALSE)
} else {
0
}
do_stop <- tox_prob_first_dose > stopping@prob
# generate message
msg <- if (n == 0) {
"Lowest active dose not tested, stopping rule not applied."
} else {
paste(
"Probability that the lowest active dose of ",
data@doseGrid[data@placebo + 1],
" being toxic based on posterior Beta distribution using a Beta(",
stopping@a, ",", stopping@b, ") prior is ",
round(tox_prob_first_dose * 100),
"% and thus ",
ifelse(do_stop, "above", "below"),
" the required ",
round(stopping@prob * 100),
"% threshold.",
sep = ""
)
}
structure(
do_stop,
message = msg,
report_label = stopping@report_label
)
}
)
## StoppingTargetBiomarker ----
#' @describeIn stopTrial Stop based on probability of targeting biomarker
#'
#' @aliases stopTrial-StoppingTargetBiomarker
#' @example examples/Rules-method-stopTrial-StoppingTargetBiomarker.R
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingTargetBiomarker",
dose = "numeric",
samples = "Samples",
model = "DualEndpoint",
data = "ANY"
),
definition = function(stopping, dose, samples, model, data, ...) {
# Compute the target biomarker prob at this dose.
# Get the biomarker level samples at the dose grid points.
biom_level_samples <- biomarker(xLevel = seq_len(data@nGrid), model, samples, ...)
# If target is relative to maximum.
if (stopping@is_relative) {
# If there is an 'Emax' parameter, target biomarker level will
# be relative to 'Emax', otherwise will be relative to the
# maximum biomarker level achieved in the given dose range.
if ("Emax" %in% names(samples)) {
# For each sample, look which dose is maximizing the
# simultaneous probability to be in the target biomarker
# range and below overdose toxicity.
prob_target <- numeric(ncol(biom_level_samples))
prob_target <- sapply(
seq(1, ncol(biom_level_samples)),
function(x) {
sum(biom_level_samples[, x] >= stopping@target[1] * samples@data$Emax &
biom_level_samples[, x] <= stopping@target[2] * samples@data$Emax) /
nrow(biom_level_samples)
}
)
} else {
# For each sample, look which was the minimum dose giving
# relative target level.
targetIndex <- apply(
biom_level_samples, 1L,
function(x) {
rnx <- range(x)
min(which(
(x >= stopping@target[1] * diff(rnx) + rnx[1]) &
(x <= stopping@target[2] * diff(rnx) + rnx[1] + 1e-10)
))
}
)
prob_target <- numeric(ncol(biom_level_samples))
tab <- table(targetIndex)
prob_target[as.numeric(names(tab))] <- tab
prob_target <- prob_target / nrow(biom_level_samples)
}
} else {
# Otherwise the target is absolute.
# For each sample, look which dose is maximizing the
# simultaneous probability to be in the target biomarker
# range and below overdose toxicity.
prob_target <- numeric(ncol(biom_level_samples))
prob_target <- sapply(
seq(1, ncol(biom_level_samples)),
function(x) {
sum(biom_level_samples[, x] >= stopping@target[1] & biom_level_samples[, x] <= stopping@target[2]) /
nrow(biom_level_samples)
}
)
}
prob_target <- ifelse(
is.na(dose),
0,
prob_target[which(data@doseGrid == dose)]
)
do_stop <- prob_target >= stopping@prob
msg <- paste(
"Probability for target biomarker is",
round(prob_target * 100),
"% for dose",
dose,
"and thus",
ifelse(do_stop, "above", "below"),
"the required",
round(stopping@prob * 100),
"%"
)
structure(
do_stop,
message = msg,
report_label = stopping@report_label
)
}
)
## StoppingSpecificDose ----
#' @describeIn stopTrial if Stopping rule is met for specific dose of the planned
#' dose grid and not just for the default next best dose.
#'
#' @aliases stopTrial-StoppingSpecificDose
#'
#' @export
#' @example examples/Rules-method-stopTrial-StoppingSpecificDose.R
#'
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingSpecificDose",
dose = "numeric",
samples = "ANY",
model = "ANY",
data = "Data"
),
definition = function(stopping, dose, samples, model, data, ...) {
# Specific dose must be a part of the dose grid.
assert_subset(x = stopping@dose@.Data, choices = data@doseGrid)
# Evaluate the original (wrapped) stopping rule at the specific dose.
result <- stopTrial(
stopping = stopping@rule,
dose = stopping@dose@.Data,
samples = samples,
model = model,
data = data,
...
)
# Correct the text message from the original stopping rule.
attr(result, "message") <- gsub(
pattern = "next best",
replacement = "specific",
x = attr(result, "message"),
ignore.case = TRUE
)
attr(result, "report_label") <- stopping@report_label
result
}
)
# nolint start
## --------------------------------------------------
## Stopping when the highest dose is reached
## --------------------------------------------------
##' @describeIn stopTrial Stop when the highest dose is reached
##'
##' @example examples/Rules-method-stopTrial-StoppingHighestDose.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingHighestDose",
dose = "numeric",
samples = "ANY",
model = "ANY",
data = "Data"
),
def =
function(stopping, dose, samples, model, data, ...) {
isHighestDose <- ifelse(
is.na(dose),
FALSE,
(dose == data@doseGrid[data@nGrid])
)
return(structure(isHighestDose,
message =
paste(
"Next best dose is", dose, "and thus",
ifelse(isHighestDose, "the",
"not the"
),
"highest dose"
),
report_label = stopping@report_label
))
}
)
## StoppingOrdinal ----
#' @describeIn stopTrial Stop based on value returned by next best dose.
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' @aliases stopTrial-StoppingOrdinal
#' @example examples/Rules-method-stopTrial-StoppingOrdinal.R
#'
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingOrdinal",
dose = "numeric",
samples = "ANY",
model = "LogisticLogNormalOrdinal",
data = "DataOrdinal"
),
definition = function(stopping, dose, samples, model, data, ...) {
stopTrial(
stopping = stopping@rule,
dose = dose,
samples = h_convert_ordinal_samples(samples, stopping@grade),
model = h_convert_ordinal_model(model, stopping@grade),
data = h_convert_ordinal_data(data, stopping@grade),
...
)
}
)
#' @describeIn stopTrial Stop based on value returned by next best dose.
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' @aliases stopTrial-StoppingOrdinal
#' @example examples/Rules-method-stopTrial-StoppingOrdinal.R
#'
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingOrdinal",
dose = "numeric",
samples = "ANY",
model = "ANY",
data = "ANY"
),
definition = function(stopping, dose, samples, model, data, ...) {
stop(
paste0(
"StoppingOrdinal objects can only be used with LogisticLogNormalOrdinal ",
"models and DataOrdinal data objects. In this case, the model is a '",
class(model),
"' object and the data is in a ",
class(data),
" object."
)
)
}
)
## StoppingExternal ----
#' @describeIn stopTrial Stop based on an external flag.
#'
#' @description `r lifecycle::badge("experimental")`
#' @param external (`flag`)\cr whether to stop based on the external
#' result or not.
#'
#' @aliases stopTrial-StoppingExternal
#' @example examples/Rules-method-stopTrial-StoppingExternal.R
#'
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingExternal",
dose = "numeric",
samples = "ANY",
model = "ANY",
data = "ANY"
),
definition = function(stopping, dose, samples, model, data, external, ...) {
assert_flag(external)
msg <- paste(
"Based on external result",
ifelse(external, "stop", "continue")
)
structure(
external,
message = msg,
report_label = stopping@report_label
)
}
)
## ============================================================
## --------------------------------------------------
## "MAX" combination of cohort size rules
## --------------------------------------------------
##' "MAX" combination of cohort size rules
##'
##' This function combines cohort size rules by taking
##' the maximum of all sizes.
##'
##' @param \dots Objects of class \code{\linkS4class{CohortSize}}
##' @return the combination as an object of class
##' \code{\linkS4class{CohortSizeMax}}
##'
##' @seealso \code{\link{minSize}}
##' @export
##' @keywords methods
setGeneric("maxSize",
def =
function(...) {
## there should be no default method,
## therefore just forward to next method!
standardGeneric("maxSize")
},
valueClass = "CohortSizeMax"
)
##' @describeIn maxSize The method combining cohort size rules by taking maximum
##' @example examples/Rules-method-maxSize.R
setMethod("maxSize",
"CohortSize",
def =
function(...) {
CohortSizeMax(list(...))
}
)
## --------------------------------------------------
## "MIN" combination of cohort size rules
## --------------------------------------------------
##' "MIN" combination of cohort size rules
##'
##' This function combines cohort size rules by taking
##' the minimum of all sizes.
##'
##' @param \dots Objects of class \code{\linkS4class{CohortSize}}
##' @return the combination as an object of class
##' \code{\linkS4class{CohortSizeMin}}
##'
##' @seealso \code{\link{maxSize}}
##' @export
##' @keywords methods
setGeneric("minSize",
def =
function(...) {
## there should be no default method,
## therefore just forward to next method!
standardGeneric("minSize")
},
valueClass = "CohortSizeMin"
)
##' @describeIn minSize The method combining cohort size rules by taking minimum
##' @example examples/Rules-method-minSize.R
setMethod("minSize",
"CohortSize",
def =
function(...) {
CohortSizeMin(list(...))
}
)
# size ----
## CohortSizeRange ----
#' @describeIn size Determines the size of the next cohort based on the range
#' into which the next dose falls into.
#'
#' @param dose the next dose.
#' @param data the data input, an object of class [`Data`].
#'
#' @aliases size-CohortSizeRange
#' @example examples/Rules-method-size-CohortSizeRange.R
#'
setMethod(
f = "size",
signature = signature(
object = "CohortSizeRange"
),
definition = function(object, dose, data) {
# If the recommended next dose is NA, don't check it and return 0.
if (is.na(dose)) {
return(0L)
}
assert_class(data, "Data")
# Determine in which interval the next dose is.
interval <- findInterval(x = dose, vec = object@intervals)
object@cohort_size[interval]
}
)
## CohortSizeDLT ----
#' @describeIn size Determines the size of the next cohort based on the number
#' of DLTs so far.
#'
#' @param dose the next dose.
#' @param data the data input, an object of class [`Data`].
#'
#' @aliases size-CohortSizeDLT
#' @example examples/Rules-method-size-CohortSizeDLT.R
#'
setMethod(
f = "size",
signature = signature(
object = "CohortSizeDLT"
),
definition = function(object, dose, data) {
# If the recommended next dose is NA, don't check it and return 0.
if (is.na(dose)) {
return(0L)
}
assert_class(data, "Data")
# Determine how many DLTs have occurred so far.
dlt_happened <- sum(data@y)
# Determine in which interval this is.
interval <- findInterval(x = dlt_happened, vec = object@intervals)
object@cohort_size[interval]
}
)
## CohortSizeMax ----
#' @describeIn size Determines the size of the next cohort based on maximum of
#' multiple cohort size rules.
#'
#' @param dose the next dose.
#' @param data the data input, an object of class [`Data`].
#'
#' @aliases size-CohortSizeMax
#' @example examples/Rules-method-size-CohortSizeMax.R
#'
setMethod(
f = "size",
signature = signature(
object = "CohortSizeMax"
),
definition = function(object, dose, data) {
# If the recommended next dose is NA, don't check it and return 0.
if (is.na(dose)) {
return(0L)
}
assert_class(data, "Data")
# Evaluate the individual cohort size rules in the list.
individual_results <- sapply(
object@cohort_sizes,
size,
dose = dose,
data = data
)
# The overall result.
max(individual_results)
}
)
## CohortSizeMin ----
#' @describeIn size Determines the size of the next cohort based on minimum of
#' multiple cohort size rules.
#'
#' @param dose the next dose.
#' @param data the data input, an object of class [`Data`].
#'
#' @aliases size-CohortSizeMin
#' @example examples/Rules-method-size-CohortSizeMin.R
#'
setMethod(
f = "size",
signature = signature(
object = "CohortSizeMin"
),
definition = function(object, dose, data) {
# If the recommended next dose is NA, don't check it and return 0.
if (is.na(dose)) {
return(0L)
}
assert_class(data, "Data")
# Evaluate the individual cohort size rules in the list.
individual_results <- sapply(
object@cohort_sizes,
size,
dose = dose,
data = data
)
# The overall result.
min(individual_results)
}
)
## CohortSizeConst ----
#' @describeIn size Constant cohort size.
#'
#' @param dose the next dose.
#' @param ... not used.
#'
#' @aliases size-CohortSizeConst
#' @example examples/Rules-method-size-CohortSizeConst.R
#'
setMethod(
f = "size",
signature = signature(
object = "CohortSizeConst"
),
definition = function(object, dose, ...) {
# If the recommended next dose is NA, don't check it and return 0.
if (is.na(dose)) {
0L
} else {
object@size
}
}
)
## CohortSizeParts ----
#' @describeIn size Determines the size of the next cohort based on the parts.
#'
#' @param dose the next dose.
#' @param data the data input, an object of class [`Data`].
#'
#' @aliases size-CohortSizeParts
#' @example examples/Rules-method-size-CohortSizeParts.R
#'
setMethod(
f = "size",
signature = signature(
object = "CohortSizeParts"
),
definition = function(object, dose, data) {
# If the recommended next dose is NA, don't check it and return 0.
if (is.na(dose)) {
return(0L)
} else {
assert_class(data, "DataParts")
object@cohort_sizes[data@nextPart]
}
}
)
## CohortSizeOrdinal ----
#' @describeIn size Determines the size of the next cohort in a ordinal CRM trial.
#'
#' @param dose (`numeric`) the next dose.
#' @param data the data input, an object of class [`DataOrdinal`].
#'
#' @aliases size-CohortSizeOrdinal
#' @example examples/Rules-method-size-CohortSizeOrdinal.R
#'
setMethod(
f = "size",
signature = signature(
object = "CohortSizeOrdinal"
),
definition = function(object, dose, data, ...) {
# Validate
assert_numeric(dose, len = 1, lower = 0)
assert_class(data, "DataOrdinal")
# Execute
size(
object@rule,
dose = dose,
data = h_convert_ordinal_data(data, object@grade),
...
)
}
)
## ------------------------------------------------------------------------------------------------
## Stopping based on a target ratio of the upper to the lower 95% credibility interval
## ------------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on 'StoppingTDCIRatio' class when
##' reaching the target ratio of the upper to the lower 95% credibility
##' interval of the estimate (TDtargetEndOfTrial). This is a stopping rule which incorporate only
##' DLE responses and DLE samples are given
##'
##' @example examples/Rules-method-stopTrialCITDsamples.R
##'
##' @export
##' @keywords methods
setMethod(
f = "stopTrial",
signature = signature(
stopping = "StoppingTDCIRatio",
dose = "ANY",
samples = "Samples",
model = "ModelTox",
data = "ANY"
),
definition = function(stopping, dose, samples, model, data, ...) {
assert_probability(stopping@prob_target)
dose_target_samples <- dose(
x = stopping@prob_target,
model = model,
samples = samples,
...
)
# 95% credibility interval.
dose_target_ci <- quantile(dose_target_samples, probs = c(0.025, 0.975))
dose_target_ci_ratio <- dose_target_ci[[2]] / dose_target_ci[[1]]
do_stop <- dose_target_ci_ratio <= stopping@target_ratio
text <- paste0(
"95% CI is (",
paste(dose_target_ci, collapse = ", "),
"), Ratio = ",
round(dose_target_ci_ratio, 4),
" is ",
ifelse(do_stop, "less than or equal to ", "greater than "),
"target_ratio = ", stopping@target_ratio
)
structure(do_stop,
message = text,
report_label = stopping@report_label
)
}
)
## ----------------------------------------------------------------------------------------------
## Stopping based on a target ratio of the upper to the lower 95% credibility interval
## ------------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on 'StoppingTDCIRatio' class
##' when reaching the target ratio of the upper to the lower 95% credibility
##' interval of the estimate (TDtargetEndOfTrial). This is a stopping rule which incorporate only
##' DLE responses and no DLE samples are involved
##' @example examples/Rules-method-stopTrialCITD.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingTDCIRatio",
dose = "ANY",
samples = "missing",
model = "ModelTox",
data = "ANY"
),
def =
function(stopping, dose, model, data, ...) {
assert_probability(stopping@prob_target)
prob_target <- stopping@prob_target
dose_target_samples <- dose(x = prob_target, model = model, ...)
## Find the variance of the log of the dose_target_samples(eta)
M1 <- matrix(c(-1 / (model@phi2), -(log(prob_target / (1 - prob_target)) - model@phi1) / (model@phi2)^2), 1, 2)
M2 <- model@Pcov
varEta <- as.vector(M1 %*% M2 %*% t(M1))
## Find the upper and lower limit of the 95% credibility interval
CI <- exp(log(dose_target_samples) + c(-1, 1) * 1.96 * sqrt(varEta))
ratio <- CI[2] / CI[1]
## so can we stop?
doStop <- ratio <= stopping@target_ratio
## generate message
text <- paste(
"95% CI is (", round(CI[1], 4), ",", round(CI[2], 4), "), Ratio =", round(ratio, 4), "is ", ifelse(doStop, "is less than or equal to", "greater than"),
"target_ratio =", stopping@target_ratio
)
## return both
return(structure(
doStop,
message = text,
report_label = stopping@report_label
))
}
)
## --------------------------------------------------------------------------------------------------
## Stopping based on a target ratio of the upper to the lower 95% credibility interval
## ------------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on reaching the target ratio of the upper to the lower 95% credibility
##' interval of the estimate (the minimum of Gstar and TDtargetEndOfTrial). This is a stopping rule which
##' incorporate DLE and efficacy responses and DLE and efficacy samples are also used.
##'
##' @param TDderive the function which derives from the input, a vector of the posterior samples called
##' \code{TDsamples} of the dose
##' which has the probability of the occurrence of DLE equals to either the targetDuringTrial or
##' targetEndOfTrial, the final next best TDtargetDuringTrial (the dose with probability of the
##' occurrence of DLE equals to the targetDuringTrial)and TDtargetEndOfTrial estimate.
##' @param Effmodel the efficacy model of \code{\linkS4class{ModelEff}} class object
##' @param Effsamples the efficacy samples of \code{\linkS4class{Samples}} class object
##' @param Gstarderive the function which derives from the input, a vector of the posterior Gstar (the dose
##' which gives the maximum gain value) samples
##' called \code{Gstarsamples}, the final next best Gstar estimate.
##'
##' @example examples/Rules-method-stopTrialCIMaxGainSamples.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingMaxGainCIRatio",
dose = "ANY",
samples = "Samples",
model = "ModelTox",
data = "DataDual"
),
def =
function(stopping, dose, samples, model, data, TDderive, Effmodel, Effsamples, Gstarderive, ...) {
prob_target <- stopping@prob_target
## checks
assert_probability(prob_target)
stopifnot(is(Effmodel, "ModelEff"))
stopifnot(is(Effsamples, "Samples"))
stopifnot(is.function(TDderive))
stopifnot(is.function(Gstarderive))
## find the TDtarget End of Trial samples
TDtargetEndOfTrialSamples <- dose(
x = prob_target,
model = model,
samples = samples,
...
)
## Find the TDtarget End of trial estimate
TDtargetEndOfTrialEstimate <- TDderive(TDtargetEndOfTrialSamples)
## Find the gain value samples then the GstarSamples
points <- data@doseGrid
GainSamples <- matrix(
nrow = size(samples),
ncol = length(points)
)
## evaluate the probs, for all gain samples.
for (i in seq_along(points))
{
## Now we want to evaluate for the
## following dose:
GainSamples[, i] <- gain(
dose = points[i],
model,
samples,
Effmodel,
Effsamples,
...
)
}
## Find the maximum gain value samples
MaxGainSamples <- apply(GainSamples, 1, max)
## Obtain Gstar samples, samples for the dose level which gives the maximum gain value
IndexG <- apply(GainSamples, 1, which.max)
GstarSamples <- data@doseGrid[IndexG]
## Find the Gstar estimate
Gstar <- Gstarderive(GstarSamples)
## Find the 95% credibility interval of Gstar and its ratio of the upper to the lower limit
CIGstar <- quantile(GstarSamples, probs = c(0.025, 0.975))
ratioGstar <- as.numeric(CIGstar[2] / CIGstar[1])
## Find the 95% credibility interval of TDtargetEndOfTrial and its ratio of the upper to the lower limit
CITDEOT <- quantile(TDtargetEndOfTrialSamples, probs = c(0.025, 0.975))
ratioTDEOT <- as.numeric(CITDEOT[2] / CITDEOT[1])
## Find which is smaller (TDtargetEndOfTrialEstimate or Gstar)
if (TDtargetEndOfTrialEstimate <= Gstar) {
## Find the upper and lower limit of the 95% credibility interval and its ratio of the smaller
CI <- CITDEOT
ratio <- ratioTDEOT
chooseTD <- TRUE
} else {
CI <- CIGstar
ratio <- ratioGstar
chooseTD <- FALSE
}
## so can we stop?
doStop <- ratio <= stopping@target_ratio
## generate message
text1 <- paste(
"Gstar estimate is", round(Gstar, 4), "with 95% CI (", round(CIGstar[1], 4), ",", round(CIGstar[2], 4),
") and its ratio =",
round(ratioGstar, 4)
)
text2 <- paste(
"TDtargetEndOfTrial estimate is ", round(TDtargetEndOfTrialEstimate, 4),
"with 95% CI (", round(CITDEOT[1], 4), ",", round(CITDEOT[2], 4), ") and its ratio=",
round(ratioTDEOT, 4)
)
text3 <- paste(
ifelse(chooseTD, "TDtargetEndOfTrial estimate", "Gstar estimate"), "is smaller with ratio =",
round(ratio, 4), " which is ", ifelse(doStop, "is less than or equal to", "greater than"),
"target_ratio =", stopping@target_ratio
)
text <- c(text1, text2, text3)
## return both
return(structure(doStop,
message = text,
report_label = stopping@report_label
))
}
)
## -----------------------------------------------------------------------------------------------
## Stopping based on a target ratio of the upper to the lower 95% credibility interval
## --------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on reaching the target ratio of the upper to the lower 95% credibility
##' interval of the estimate (the minimum of Gstar and TDtargetEndOfTrial). This is a stopping rule which
##' incorporate DLE and efficacy responses without DLE and efficacy samples involved.
##' @example examples/Rules-method-stopTrialCIMaxGain.R
setMethod("stopTrial",
signature =
signature(
stopping = "StoppingMaxGainCIRatio",
dose = "ANY",
samples = "missing",
model = "ModelTox",
data = "DataDual"
),
def =
function(stopping, dose, model, data, Effmodel, ...) {
prob_target <- stopping@prob_target
## checks
assert_probability(prob_target)
stopifnot(is(Effmodel, "ModelEff"))
## find the TDtarget End of Trial
TDtargetEndOfTrial <- dose(
x = prob_target,
model = model,
...
)
## Find the dose with maximum gain value
Gainfun <- function(DOSE) {
-gain(DOSE, model_dle = model, model_eff = Effmodel, ...)
}
# if(data@placebo) {
# n <- length(data@doseGrid)
# LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
# }
Gstar <- (optim(LowestDose, Gainfun, method = "L-BFGS-B", lower = LowestDose, upper = max(data@doseGrid))$par)
MaxGain <- -(optim(LowestDose, Gainfun, method = "L-BFGS-B", lower = LowestDose, upper = max(data@doseGrid))$value)
if (data@placebo) {
logGstar <- log(Gstar + Effmodel@const)
} else {
logGstar <- log(Gstar)
}
## From paper (Yeung et. al 2015)
meanEffGstar <- Effmodel@theta1 + Effmodel@theta2 * log(logGstar)
denom <- (model@phi2) * (meanEffGstar) * (1 + logGstar * model@phi2)
dgphi1 <- -(meanEffGstar * logGstar * model@phi2 - Effmodel@theta2) / denom
dgphi2 <- -((meanEffGstar) * logGstar + meanEffGstar * (logGstar)^2 * model@phi2 - Effmodel@theta2 * logGstar) / denom
dgtheta1 <- -(logGstar * model@phi2) / denom
dgtheta2 <- -(logGstar * exp(model@phi1 + model@phi2 * logGstar) * model@phi2 * log(logGstar) - 1 - exp(model@phi1 + model@phi2 * logGstar)) / denom
# DLEPRO <- exp(model@phi1+model@phi2*logGstar)
# dgphi1 <- Effmodel@theta2*DLEPRO - logGstar*model@phi2*meanEffGstar*DLEPRO
# dgphi2 <- logGstar*DLEPRO *(Effmodel@theta2-(meanEffGstar)+model@phi2)
# dgtheta1 <- -logGstar*DLEPRO*model@phi2
# dgtheta2 <- 1+DLEPRO-logGstar*DLEPRO*model@phi2*log(logGstar)
deltaG <- matrix(c(dgphi1, dgphi2, dgtheta1, dgtheta2), 4, 1)
## Find the variance of the log Gstar
## First find the covariance matrix of all the parameters, phi1, phi2, theta1 and theta2
## such that phi1 and phi2 and independent of theta1 and theta2
emptyMatrix <- matrix(0, 2, 2)
covBETA <- cbind(rbind(model@Pcov, emptyMatrix), rbind(emptyMatrix, Effmodel@Pcov))
varlogGstar <- as.vector(t(deltaG) %*% covBETA %*% deltaG)
## Find the upper and lower limit of the 95% credibility interval of Gstar
CIGstar <- exp(logGstar + c(-1, 1) * 1.96 * sqrt(varlogGstar))
## The ratio of the upper to the lower 95% credibility interval
ratioGstar <- CIGstar[2] / CIGstar[1]
## Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1 / (model@phi2), -(log(prob_target / (1 - prob_target)) - model@phi1) / (model@phi2)^2), 1, 2)
M2 <- model@Pcov
varEta <- as.vector(M1 %*% M2 %*% t(M1))
## Find the upper and lower limit of the 95% credibility interval of
## TDtargetEndOfTrial
CITDEOT <- exp(log(TDtargetEndOfTrial) + c(-1, 1) * 1.96 * sqrt(varEta))
## The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- CITDEOT[2] / CITDEOT[1]
if (Gstar <= TDtargetEndOfTrial) {
chooseTD <- FALSE
CI <- c()
CI[2] <- CIGstar[2]
CI[1] <- CIGstar[1]
ratio <- ratioGstar
} else {
chooseTD <- TRUE
CI <- c()
CI[2] <- CITDEOT[2]
CI[1] <- CITDEOT[1]
ratio <- ratioTDEOT
}
## so can we stop?
doStop <- ratio <= stopping@target_ratio
## generate message
text1 <- paste(
"Gstar estimate is", round(Gstar, 4), "with 95% CI (", round(CIGstar[1], 4), ",", round(CIGstar[2], 4),
") and its ratio =",
round(ratioGstar, 4)
)
text2 <- paste(
"TDtargetEndOfTrial estimate is ", round(TDtargetEndOfTrial, 4),
"with 95% CI (", round(CITDEOT[1], 4), ",", round(CITDEOT[2], 4), ") and its ratio=",
round(ratioTDEOT, 4)
)
text3 <- paste(
ifelse(chooseTD, "TDatrgetEndOfTrial estimate", "Gstar estimate"), "is smaller with ratio =",
round(ratio, 4), "which is ", ifelse(doStop, "is less than or equal to", "greater than"),
"target_ratio =", stopping@target_ratio
)
text <- c(text1, text2, text3)
## return both
return(structure(doStop,
message = text,
report_label = stopping@report_label
))
}
)
## ============================================================
## -----------------------------------------------------
## Determine the safety window length of the next cohort
## -----------------------------------------------------
##' Determine the safety window length of the next cohort
##'
##' This function determines the safety window length of
##' the next cohort.
##'
##' @param safetyWindow The rule, an object of class
##' \code{\linkS4class{SafetyWindow}}
##' @param size The next cohort size
##' @param data The data input, an object of class \code{\linkS4class{DataDA}}
##' @param \dots additional arguments
##'
##' @return the `windowLength` as a list of safety window parameters
##' (`gap`, `follow`, `follow_min`)
##'
##' @export
##' @keywords methods
setGeneric("windowLength",
def =
function(safetyWindow, size, ...) {
## there should be no default method,
## therefore just forward to next method!
standardGeneric("windowLength")
},
valueClass = "list"
)
## ============================================================
## --------------------------------------------------
## The SafetyWindowSize method
## --------------------------------------------------
##' @describeIn windowLength Determine safety window length based
##' on the cohort size
##'
##' @example examples/Rules-method-windowLength-SafetyWindowSize.R
setMethod("windowLength",
signature =
signature(
safetyWindow = "SafetyWindowSize",
size = "ANY"
),
def =
function(safetyWindow, size, data, ...) {
## determine in which interval the next size is
interval <-
findInterval(
x = size,
vec = safetyWindow@size
)
## so the safety window length is
patientGap <- head(c(
0, safetyWindow@gap[[interval]],
rep(tail(safetyWindow@gap[[interval]], 1), 100)
), size)
patientFollow <- safetyWindow@follow
patientFollowMin <- safetyWindow@follow_min
ret <- list(patientGap = patientGap, patientFollow = patientFollow, patientFollowMin = patientFollowMin)
return(ret)
}
)
## ============================================================
## --------------------------------------------------
## Constant safety window length
## --------------------------------------------------
##' @describeIn windowLength Constant safety window length
##' @example examples/Rules-method-windowLength-SafetyWindowConst.R
setMethod("windowLength",
signature =
signature(
safetyWindow = "SafetyWindowConst",
size = "ANY"
),
def =
function(safetyWindow, size, ...) {
## first element should be 0.
patientGap <- head(c(
0, safetyWindow@gap,
rep(tail(safetyWindow@gap, 1), 100)
), size)
patientFollow <- safetyWindow@follow
patientFollowMin <- safetyWindow@follow_min
ret <- list(
patientGap = patientGap,
patientFollow = patientFollow,
patientFollowMin = patientFollowMin
)
return(ret)
}
)
# nolint end
# tidy ----
## tidy-IncrementsRelative ----
#' @rdname tidy
#' @aliases tidy-IncrementsRelative
#' @example examples/Rules-method-tidyIncrementsRelative.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "IncrementsRelative"),
definition = function(x, ...) {
h_tidy_all_slots(x) %>%
dplyr::bind_cols() %>%
h_range_to_minmax(.data$intervals) %>%
dplyr::filter(max > 0) %>%
tibble::add_column(increment = x@increments) %>%
h_tidy_class(x)
}
)
## tidy-CohortSizeDLT ----
#' @rdname tidy
#' @aliases tidy-CohortSizeDLT
#' @example examples/Rules-method-tidyCohortSizeDLT.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "CohortSizeDLT"),
definition = function(x, ...) {
h_tidy_all_slots(x) %>%
dplyr::bind_cols() %>%
h_range_to_minmax(.data$intervals) %>%
dplyr::filter(max > 0) %>%
tibble::add_column(cohort_size = x@cohort_size) %>%
h_tidy_class(x)
}
)
## tidy-CohortSizeMin ----
#' @rdname tidy
#' @aliases tidy-CohortSizeMin
#' @example examples/Rules-method-tidyCohortSizeMin.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "CohortSizeMin"),
definition = function(x, ...) {
callNextMethod() %>% h_tidy_class(x)
}
)
## tidy-CohortSizeMax ----
#' @rdname tidy
#' @aliases tidy-CohortSizeMax
#' @example examples/Rules-method-tidyCohortSizeMax.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "CohortSizeMax"),
definition = function(x, ...) {
callNextMethod() %>% h_tidy_class(x)
}
)
## tidy-CohortSizeRange ----
#' @rdname tidy
#' @aliases tidy-CohortSizeRange
#' @example examples/Rules-method-tidyCohortSizeRange.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "CohortSizeRange"),
definition = function(x, ...) {
h_tidy_all_slots(x) %>%
dplyr::bind_cols() %>%
h_range_to_minmax(.data$intervals) %>%
dplyr::filter(max > 0) %>%
tibble::add_column(cohort_size = x@cohort_size) %>%
h_tidy_class(x)
}
)
## tidy-CohortSizeParts ----
#' @rdname tidy
#' @aliases tidy-CohortSizeParts
#' @example examples/Rules-method-tidyCohortSizeParts.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "CohortSizeParts"),
definition = function(x, ...) {
tibble::tibble(
part = seq_along(x@cohort_sizes),
cohort_size = x@cohort_sizes
) %>%
h_tidy_class(x)
}
)
## tidy-IncrementsMin ----
#' @rdname tidy
#' @aliases tidy-IncrementsMin
#' @example examples/Rules-method-tidyIncrementsMin.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "IncrementsMin"),
definition = function(x, ...) {
callNextMethod() %>% h_tidy_class(x)
}
)
## tidy-IncrementsRelative ----
#' @rdname tidy
#' @aliases tidy-IncrementsRelative
#' @example examples/Rules-method-tidyIncrementsRelative.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "IncrementsRelative"),
definition = function(x, ...) {
h_tidy_all_slots(x) %>%
h_range_to_minmax(.data$intervals) %>%
dplyr::filter(dplyr::row_number() > 1) %>%
tibble::add_column(increment = x@increments) %>%
h_tidy_class(x)
}
)
## tidy-IncrementsRelativeDLT ----
#' @rdname tidy
#' @aliases tidy-IncrementsRelativeDLT
#' @example examples/Rules-method-tidyIncrementsRelativeDLT.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "IncrementsRelativeDLT"),
definition = function(x, ...) {
h_tidy_all_slots(x) %>%
h_range_to_minmax(.data$intervals) %>%
dplyr::filter(dplyr::row_number() > 1) %>%
tibble::add_column(increment = x@increments) %>%
h_tidy_class(x)
}
)
## tidy-IncrementsRelative ----
#' @rdname tidy
#' @aliases tidy-IncrementsRelativeParts
#' @example examples/Rules-method-tidyIncrementsRelativeParts.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "IncrementsRelativeParts"),
definition = function(x, ...) {
slot_names <- slotNames(x)
rv <- list()
for (nm in slot_names) {
if (!is.function(slot(x, nm))) {
rv[[nm]] <- h_tidy_slot(x, nm, ...)
}
}
# Column bind of all list elements have the same number of rows.
if (length(rv) > 1 & length(unique(sapply(rv, nrow))) == 1) {
rv <- rv %>% dplyr::bind_cols()
}
rv <- rv %>% h_tidy_class(x)
if (length(rv) == 1) {
rv[[names(rv)[1]]] %>% h_tidy_class(x)
} else {
rv
}
}
)
## tidy-NextBestNCRM ----
#' @rdname tidy
#' @aliases tidy-NextBestNCRM
#' @example examples/Rules-method-tidyNextBestNCRM.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "NextBestNCRM"),
definition = function(x, ...) {
h_tidy_all_slots(x) %>%
dplyr::bind_cols() %>%
h_range_to_minmax(.data$target, range_min = 0, range_max = 1) %>%
add_column(max_prob = c(NA, NA, x@max_overdose_prob)) %>%
add_column(Range = c("Underdose", "Target", "Overdose"), .before = 1) %>%
h_tidy_class(x)
}
)
## tidy-NextBestNCRMLoss ----
#' @rdname tidy
#' @aliases tidy-NextBestNCRMLoss
#' @example examples/Rules-method-tidyNextBestNCRMLoss.R
#' @export
setMethod(
f = "tidy",
signature = signature(x = "NextBestNCRMLoss"),
definition = function(x, ...) {
tibble(
Range = "Underdose",
Lower = 0,
Upper = x@target[1]
) %>%
dplyr::bind_rows(
lapply(
c("target", "overdose", "unacceptable"),
function(nm, obj) {
tibble::tibble(
Range = stringr::str_to_sentence(nm),
Lower = slot(obj, nm)[1],
Upper = slot(obj, nm)[2]
)
},
obj = x
) %>% dplyr::bind_rows()
) %>%
add_column(LossCoefficient = x@losses) %>%
add_column(MaxOverdoseProb = x@max_overdose_prob) %>%
h_tidy_class(x)
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.