Nothing
R/get.oc.comb.R
In BOIN: Bayesian Optimal INterval (BOIN) Design for Single-Agent and Drug- Combination Phase I Clinical Trials
#'
#' Generate operating characteristics for drug combination trials
#'
#' Obtain the operating characteristics of the BOIN design or waterfall design for drug combination
#' trials. The BOIN design is to find a MTD, and the waterfall design is to find the MTD contour
#' (i.e., multple MTDs in the dose matrix)
#'
#' @usage get.oc.comb(target, p.true, ncohort, cohortsize, n.earlystop=NULL, startdose=c(1, 1),
#' titration=FALSE,p.saf=0.6*target, p.tox=1.4*target, cutoff.eli=0.95,
#' extrasafe=FALSE,offset=0.05, ntrial=1000, mtd.contour=FALSE,
#' boundMTD=FALSE, seed=6)
#'
#' @param target the target DLT rate
#' @param p.true a \code{J*K} matrix \code{(J<=K)} containing the true toxicity probabilities of
#' combinations with \code{J} dose levels of agent A and \code{K} dose levels of agent B
#' @param ncohort a \code{1*J} vector specifying the number of cohorts for each of \code{J} subtrials
#' if \code{mtd.contour=TRUE}; Otherwise, a scalar specifying the total number of cohorts for
#' the trial.
#' @param cohortsize the cohort size
#' @param n.earlystop the early stopping parameter. If the number of patients treated at the current
#' dose reaches \code{n.earlystop}, stop the trial or subtrial and select the MTD based on
#' the observed data. When the waterfall design is used to find the MTD contour,
#' \code{n.earlystop=12} by default.
#' @param startdose the starting dose combination level for drug combination trial
#' @param titration set \code{titration=TRUE} to perform dose escalation with cohort size = 1 to accelerate dose escalation at the begining of the trial.
#' @param p.saf the highest toxicity probability that is deemed subtherapeutic (i.e. below the MTD)
#' such that dose escalation should be undertaken.
#' The default value is \code{p.saf=0.6*target}.
#' @param p.tox the lowest toxicity probability that is deemed overly toxic such that deescalation
#' is required. The default value is \code{p.tox=1.4*target}.
#' @param cutoff.eli the cutoff to eliminate an overly toxic dose for safety. We recommend the
#' default value of (\code{cutoff.eli=0.95}) for general use.
#' @param extrasafe set \code{extrasafe=TRUE} to impose a more stringent stopping rule
#' @param offset a small positive number (between 0 and 0.5) to control how strict the stopping
#' rule is when \code{extrasafe=TRUE}. A larger value leads to a more strict stopping
#' rule. The default value \code{offset=0.05} generally works well.
#' @param ntrial the total number of trials to be simulated
#' @param mtd.contour set \code{mtd.contour=TRUE} to select the MTD contour (claiming multiple MTDs).
#' Otherwise, BOIN design is used to search for a single MTD.
#' @param boundMTD set \code{boundMTD=TRUE} to impose the condition: the isotonic estimate of toxicity
#' probability for the selected MTD must be less than de-escalation boundary.
#' @param seed the random seed for simulation
#'
#' @details The operating characteristics of the BOIN design or waterfall design are generated by
#' simulating trials under the prespecified true toxicity probabilities of the investigational dose
#' combinations. If \code{titration=TRUE}, we perform dose escalation with cohort size = 1 at the begining of the trial:
#' starting from \code{startdose}, if no toxicity is observed, we escalate the dose;
#' otherwise, the titration is completed and we switch to cohort size = \code{cohortsize}.
#' Titration accelerates the dose escalation and is useful when low doses are believed to be safe.
#'
#'
#' The BOIN and waterfall designs have two built-in stopping rules:
#' (1) stop the trial/subtrial if the lowest dose is eliminated due to toxicity, and no dose should
#' be selected as the MTD; and (2) stop the trial/subtrial and select the MTD if the number of
#' patients treated at the current dose reaches \code{n.earlystop}. The first stopping rule is a safety
#' rule to protect patients from the case in which all doses are overly toxic. The rationale for
#' the second stopping rule is that when there is a large number (i.e., \code{n.earlystop}) of
#' patients assigned to a dose, it means that the dose-finding algorithm has approximately converged.
#' Thus, we can stop the trial/subtrial early and select the MTD to save sample size and reduce the
#' trial duration.
#'
#'
#' For some applications, investigators may prefer a more strict safety stopping rule than rule
#' (1) for extra safety when the lowest dose is overly toxic.
#' This can be achieved by setting \code{extrasafe=TRUE},
#' which imposes the following more strict safety stopping rule:
#' stop the trial if (i) the number of patients treated at the lowest dose \code{>=3},
#' and (ii) \eqn{Pr(toxicity\ rate\ of\ the\ lowest\ dose > \code{target} | data) > \code{cutoff.eli}-\code{offset}}.
#' As a tradeoff, the strong stopping rule will decrease the MTD selection percentage
#' when the lowest dose actually is the MTD.
#'
#' @return \code{get.oc.comb()} returns the operating characteristics of the BOIN combination or
#' waterfall design as a list. For the BOIN combination design, including:
#' (1) true toxicity probability at each dose level (\code{$p.true}),
#' (2) selection percentage at each dose level (\code{$selpercent}),
#' (3) the number of patients treated at each dose level (\code{$npatients})
#' (4) the number of toxicities observed at each dose level (\code{$ntox}),
#' (5) the total number of patients in the trial (\code{$totaln}),
#' (6) the total number of toxicities observed for the trial (\code{$totaltox})
#' (7) the pecentage of correct selection (\code{$pcs}),
#' (8) the total percentage of patients treated at the MTD (\code{$npercent}).
#' (9) the percentage of early stopping without selecting the MTD (\code{$percentstop})
#' For the the waterfall design, including:
#' (1) true toxicity probability at each dose level (\code{$p.true}),
#' (2) selection percentage of dose combinations (\code{$selpercent}),
#' (3) the number of patients treated at each dose combination (\code{$npatients})
#' (4) the number of toxicities observed at each dose combination (\code{$ntox}),
#' (5) the total number of patients in the trial (\code{$totaln}),
#' (6) the total number of toxicities observed for the trial (\code{$totaltox})
#' (7) the total percentage of correct selection at the MTD contour (\code{$pcs.contour}),
#' (8) the total percentage of patients treated at MTD contour
#' (\code{$npercent.contour})
#' (9) the total percentage of patients treated above MTD contour
#' (\code{$npercent.above.contour})
#' (10) the total percentage of patients treated below MTD contour
#' (\code{$npercent.below.contour})
#'
#'
#' @note We should avoid setting the values of \code{p.saf} and \code{p.tox} very close to the
#' \code{target}. This is because the small sample sizes of typical phase I trials prevent us from
#' differentiating the target DLT rate from the rates close to it. The default values provided by
#' \code{get.oc()} are strongly recommended, and generally yield excellent operating characteristics.
#'
#'
#' @author Suyu Liu, Liangcai Zhang, Yanhong Zhou, and Ying Yuan
#'
#' @references Liu S. and Yuan, Y. (2015). Bayesian Optimal Interval Designs for Phase I Clinical
#' Trials, \emph{Journal of the Royal Statistical Society: Series C}, 64, 507-523.
#'
#' Lin R. and Yin, G. (2017). Bayesian Optimal Interval Designs for Dose Finding in
#' Drug-combination Trials, \emph{Statistical Methods in Medical Research}, 26, 2155-2167.
#'
#' Yan, F., Zhang, L., Zhou, Y., Pan, H., Liu, S. and Yuan, Y. (2020).BOIN: An R Package
#' for Designing Single-Agent and Drug-Combination Dose-Finding Trials Using Bayesian Optimal
#' Interval Designs. \emph{Journal of Statistical Software}, 94(13),1-32.<doi:10.18637/jss.v094.i13>.
#'
#'
#' Zhang L. and Yuan, Y. (2016). A Simple Bayesian Design to Identify the Maximum
#' Tolerated Dose Contour for Drug Combination Trials, \emph{Statistics in Medicine}, 35, 4924-4936.
#'
#' @seealso Tutorial: \url{http://odin.mdacc.tmc.edu/~yyuan/Software/BOIN/BOIN2.6_tutorial.pdf}
#'
#' Paper: \url{http://odin.mdacc.tmc.edu/~yyuan/Software/BOIN/paper.pdf}
#'
#' @examples
#'
#' ###### drug-combination trial ######
#'
#' ##### combination trial to find a single MTD ######
#'
#' ## get the operating characteristics for BOIN design
#' p.true <- matrix(c(0.01,0.03,0.10,0.20,0.30,
#' 0.03,0.05,0.15,0.30,0.60,
#' 0.08,0.10,0.30,0.60,0.75), byrow=TRUE, ncol=5)
#'
#' oc.comb <- get.oc.comb(target=0.3, p.true, ncohort=20, cohortsize=3,
#' n.earlystop=12, startdose=c(1,1), ntrial=100)
#'
#' summary(oc.comb)
#' plot(oc.comb)
#'
#'
#' ## get the operating characteristics with titration for BOIN design
#' oc.comb <- get.oc.comb(target=0.3, p.true, ncohort=20, cohortsize=3,
#' n.earlystop=12, startdose=c(1,1), titration=TRUE, ntrial=100)
#' summary(oc.comb)
#' plot(oc.comb)
#'
#'
#' ##### combination trial to find the MTD contour ######
#'
#' ## find the MTD contour using waterfall design
#' oc.comb <- get.oc.comb(target=0.3, p.true, ncohort=c(10,5,5), cohortsize=3,
#' n.earlystop=12, startdose=c(1,1), ntrial=100, mtd.contour=TRUE)
#'
#' summary(oc.comb)
#' plot(oc.comb)
#'
#' @export
get.oc.comb <- function (target, p.true, ncohort, cohortsize, n.earlystop = NULL,
startdose = c(1, 1), titration = FALSE, p.saf = 0.6 * target,
p.tox = 1.4 * target, cutoff.eli = 0.95, extrasafe = FALSE,
offset = 0.05, ntrial = 1000, mtd.contour = FALSE, boundMTD=FALSE,seed = 6)
{
get.oc.comb.boin <- function(target, p.true, ncohort, cohortsize,
n.earlystop = 100, startdose = c(1, 1), titration = FALSE,
p.saf = 0.6 * target, p.tox = 1.4 * target, cutoff.eli = 0.95,
extrasafe = FALSE, offset = 0.05, boundMTD=FALSE,ntrial = 1000) {
JJ = nrow(p.true)
KK = ncol(p.true)
if (JJ > KK)
stop("p.true should be arranged in a way (i.e., rotated) such that\n the number of rows is less than or equal to the number of columns.")
if (JJ > KK)
p.true = t(p.true)
if (target < 0.05) {
stop("the target is too low!")
}
if (target > 0.6) {
stop("the target is too high!")
}
if ((target - p.saf) < (0.1 * target)) {
stop("the probability deemed safe cannot be higher than or too close to the target!")
}
if ((p.tox - target) < (0.1 * target)) {
stop("the probability deemed toxic cannot be lower than or too close to the target!")
}
if (offset >= 0.5) {
stop("the offset is too large!")
}
if (n.earlystop <= 6) {
warning("the value of n.earlystop is too low to ensure good operating characteristics. Recommend n.earlystop = 9 to 18")
}
ndose = length(p.true)
npts = ncohort * cohortsize
Y <- array(matrix(rep(0, length(p.true) * ntrial), dim(p.true)[1]),
dim = c(dim(p.true), ntrial))
N <- array(matrix(rep(0, length(p.true) * ntrial), dim(p.true)[1]),
dim = c(dim(p.true), ntrial))
dselect = matrix(rep(0, 2 * ntrial), ncol = 2)
if (cohortsize > 1) {
temp = get.boundary(target, ncohort, cohortsize,
n.earlystop=100, p.saf, p.tox, cutoff.eli, extrasafe)$full_boundary_tab
}else {
temp = get.boundary(target, ncohort, cohortsize,
n.earlystop=100, p.saf, p.tox, cutoff.eli, extrasafe)$boundary_tab
}
b.e = temp[2, ]
b.d = temp[3, ]
b.elim = temp[4, ]
lambda1 = log((1 - p.saf)/(1 - target))/log(target *
(1 - p.saf)/(p.saf * (1 - target)))
lambda2 = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -
target)/(target * (1 - p.tox)))
if (cohortsize == 1)
titration = FALSE
for (trial in 1:ntrial) {
y <- matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
n <- matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
earlystop = 0
d = startdose
elimi = matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
ft=TRUE #flag used to determine whether or not to add cohortsize-1 patients to a dose for the first time when titration is triggered.
if (titration) {
tmpa = d[1]
tmpb = d[2]
y[tmpa, tmpb] <- (runif(1) < p.true[tmpa, tmpb])
n[tmpa, tmpb] <- 1
while (tmpa <= dim(p.true)[1] & tmpb <= dim(p.true)[2]) {
if (tmpa == dim(p.true)[1] & tmpb == dim(p.true)[2]) {
break
}
if (sum(y) == 1) {
y[tmpa, tmpb] = 1
break
}
if (tmpa < dim(p.true)[1] & tmpb < dim(p.true)[2]) {
tmp.candidate = rbind(c(tmpa + 1, tmpb),
c(tmpa, tmpb + 1))
tmp.sel = rbinom(1, 1, prob = c(0.5, 0.5)) +
1
tmpa = tmp.candidate[tmp.sel, 1]
tmpb = tmp.candidate[tmp.sel, 2]
}
else if (tmpa == dim(p.true)[1]) {
tmpb = tmpb + 1
}
else {
tmpa = tmpa + 1
}
y[tmpa, tmpb] <- (runif(1) < p.true[tmpa, tmpb])
n[tmpa, tmpb] <- 1
}
if (sum(y) == 0) {
d = c(dim(p.true)[1], dim(p.true)[2])
}
else {
d = c(tmpa, tmpb)
}
}
for (pp in 1:ncohort) {
if (titration & n[d[1], d[2]] < cohortsize & ft) {
ft=FALSE
y[d[1], d[2]] = y[d[1], d[2]] + sum(runif(cohortsize -
1) < p.true[d[1], d[2]])
n[d[1], d[2]] = n[d[1], d[2]] + cohortsize -1
}
else {
y[d[1], d[2]] = y[d[1], d[2]] + sum(runif(cohortsize) <
p.true[d[1], d[2]])
n[d[1], d[2]] = n[d[1], d[2]] + cohortsize
}
nc = n[d[1], d[2]]
if (!is.na(b.elim[nc])) {
if (y[d[1], d[2]] >= b.elim[nc]) {
for (i in min(d[1], dim(p.true)[1]):dim(p.true)[1]) {
for (j in min(d[2], dim(p.true)[2]):dim(p.true)[2]) {
elimi[i, j] = 1
}
}
if (d[1] == 1 && d[2] == 1) {
d = c(99, 99)
earlystop = 1
break
}
}
if (extrasafe) {
if (d[1] == 1 && d[2] == 1 && n[1, 1] >=
3) {
if (1 - pbeta(target, y[1, 1] + 1, n[1,
1] - y[1, 1] + 1) > cutoff.eli - offset) {
d = c(99, 99)
earlystop = 1
break
}
}
}
}
if (n[d[1],d[2]] >= n.earlystop && (y[d[1],d[2]]>b.e[n[d[1],d[2]]] ||
(d[1]==dim(p.true)[1] && d[2]==dim(p.true)[2]) ||
( d[1]==dim(p.true)[1] && d[2]<dim(p.true)[2] && elimi[d[1],d[2]+1]==1 ) ||
( d[1]<dim(p.true)[1] && d[2]==dim(p.true)[2] && elimi[d[1]+1,d[2]]==1 ) ||
( d[1]<dim(p.true)[1] && d[2]<dim(p.true)[2] && elimi[d[1]+1,d[2]]==1 && elimi[d[1],d[2]+1]==1 ) ) &&
(y[d[1],d[2]]<b.d[n[d[1],d[2]]] || (d[1]==1 && d[2]==1) ) ) break;
if (y[d[1], d[2]] <= b.e[nc]) {
elevel = matrix(c(1, 0, 0, 1), 2)
pr_H0 = rep(0, length(elevel)/2)
nn = pr_H0
for (i in seq(1, length(elevel)/2, by = 1)) {
if (d[1] + elevel[1, i] <= dim(p.true)[1] &&
d[2] + elevel[2, i] <= dim(p.true)[2]) {
if (elimi[d[1] + elevel[1, i], d[2] + elevel[2,
i]] == 0) {
yn = y[d[1] + elevel[1, i], d[2] + elevel[2,
i]]
nn[i] = n[d[1] + elevel[1, i], d[2] +
elevel[2, i]]
pr_H0[i] <- pbeta(lambda2, yn + 0.5,
nn[i] - yn + 0.5) - pbeta(lambda1,
yn + 0.5, nn[i] - yn + 0.5)
}
}
}
pr_H0 = pr_H0 + nn * 5e-04
if (max(pr_H0) == 0) {
d = d
}
else {
k = which(pr_H0 == max(pr_H0))[as.integer(runif(1) *
length(which(pr_H0 == max(pr_H0))) + 1)]
d = d + c(elevel[1, k], elevel[2, k])
}
}
else if (y[d[1], d[2]] >= b.d[nc]) {
delevel = matrix(c(-1, 0, 0, -1), 2)
pr_H0 = rep(0, length(delevel)/2)
nn = pr_H0
for (i in seq(1, length(delevel)/2, by = 1)) {
if (d[1] + delevel[1, i] > 0 && d[2] + delevel[2,
i] > 0) {
yn = y[d[1] + delevel[1, i], d[2] + delevel[2,
i]]
nn[i] = n[d[1] + delevel[1, i], d[2] +
delevel[2, i]]
pr_H0[i] = pbeta(lambda2, yn + 0.5, nn[i] -
yn + 0.5) - pbeta(lambda1, yn + 0.5,
nn[i] - yn + 0.5)
}
}
pr_H0 = pr_H0 + nn * 5e-04
if (max(pr_H0) == 0) {
d = d
}
else {
k = which(pr_H0 == max(pr_H0))[as.integer(runif(1) *
length(which(pr_H0 == max(pr_H0))) + 1)]
d = d + c(delevel[1, k], delevel[2, k])
}
}
else {
d = d
}
}
Y[, , trial] = y
N[, , trial] = n
if (earlystop == 1) {
dselect[trial, ] = c(99, 99)
}else {
selcomb = select.mtd.comb.boin(target, n, y,
cutoff.eli, extrasafe, offset,
boundMTD=boundMTD,p.tox=p.tox,mtd.contour = FALSE)$MTD
dselect[trial, 1] = selcomb[1]
dselect[trial, 2] = selcomb[2]
}
}
selpercent = matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
nptsdose = apply(N, c(1, 2), mean, digits = 2, format = "f")
ntoxdose = apply(Y, c(1, 2), mean, digits = 2, format = "f")
for (i in 1:dim(p.true)[1]) for (j in 1:dim(p.true)[2]) {
{
selpercent[i, j] = sum(dselect[, 1] == i & dselect[,
2] == j)/ntrial * 100
}
}
if (JJ <= KK) {
rownames(p.true) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(p.true) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
out = list(p.true = round(p.true, 2), selpercent = round(selpercent,2),
npatients = round(apply(N, c(1, 2), mean),2), ntox = round(apply(Y, c(1, 2), mean), 2),
totaltox = round(sum(Y)/ntrial, 1), totaln = round(sum(N)/ntrial,1),
pcs = paste(round(sum(selpercent[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)]),1), "%", sep = ""),
npercent = paste(round(sum(nptsdose[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)])/sum(nptsdose) *100, 1), "%", sep = ""),
percentstop=100-sum(round(selpercent,2)),flowchart = FALSE)
rownames(out$npatients) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(out$npatients) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(out$ntox) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(out$ntox) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
return(out)
}
else {
colnames(p.true) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(p.true) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(npatients) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(npatients) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(ntox) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(ntox) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
out = list(p.true = round(t(p.true), 2), selpercent = round(t(selpercent),2),
npatients = round(t(apply(N, c(1, 2), mean)), 2), ntox = round(t(apply(Y, c(1, 2), mean)),2),
totaltox = round(sum(Y)/ntrial, 1), totaln = round(sum(N)/ntrial,1), pcs = paste(round(sum(selpercent[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)]),
1), "%"), npercent = paste(round(sum(nptsdose[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)])/sum(nptsdose) *100, 1), "%"),
percentstop=100-sum(round(selpercent,2)),
flowchart = FALSE)
return(out)
}
}
select.mtd.comb.boin <- function(target, npts, ntox, cutoff.eli = 0.95,
extrasafe = FALSE, offset = 0.05,
boundMTD=FALSE,p.tox=1.4*target,
mtd.contour = FALSE) {
lambda_d = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -target)/(target * (1 - p.tox)))
y = ntox
n = npts
if (nrow(n) > ncol(n) | nrow(y) > ncol(y)) {
stop("npts and ntox should be arranged in a way (i.e., rotated) such that for each of them, the number of rows is less than or equal to the number of columns.")
}
elimi = matrix(0, dim(n)[1], dim(n)[2])
if (extrasafe) {
if (n[1, 1] >= 3) {
if (1 - pbeta(target, y[1, 1] + 1, n[1, 1] -
y[1, 1] + 1) > cutoff.eli - offset) {
elimi[, ] = 1
}
}
}
for (i in 1:dim(n)[1]) {
for (j in 1:dim(n)[2]) {
if (n[i, j] >= 3) {
if (1 - pbeta(target, y[i, j] + 1, n[i, j] -
y[i, j] + 1) > cutoff.eli) {
elimi[i:dim(n)[1], j] = 1
elimi[i, j:dim(n)[2]] = 1
break
}
}
}
}
selectdose=NULL
if (elimi[1] == 1) {
selectdose = c(99, 99)
selectdoses = matrix(selectdose, nrow = 1)
}else {
phat = (y + 0.05)/(n + 0.1)
phat = Iso::biviso(phat, n + 0.1, warn = TRUE)[,
]
phat.out = phat
phat.out[n == 0] = NA
phat[elimi == 1] = 1.1
phat = phat * (n != 0) + (1e-05) * (matrix(rep(1:dim(n)[1],
each = dim(n)[2], len = length(n)), dim(n)[1],
byrow = T) + matrix(rep(1:dim(n)[2], each = dim(n)[1],
len = length(n)), dim(n)[1]))
if(boundMTD){
if(all(phat[n!=0]>=lambda_d)){
selectdose = c(99, 99)
selectdoses = matrix(selectdose, nrow = 1)
}else{
phat[phat>=lambda_d]=10}}
if(is.null(selectdose)){
phat[n == 0] = 10
selectdose = which(abs(phat - target) == min(abs(phat -
target)), arr.ind = TRUE)
if (length(selectdose) > 2)
selectdose = selectdose[1, ]
aa = function(x) as.numeric(as.character(x))
if (mtd.contour == TRUE) {
selectdoses = cbind(row = 1:dim(n)[1], col = rep(99,
dim(n)[1]))
for (k in dim(n)[1]:1) {
kn = n[k, ]
ky = y[k, ]
kelimi = elimi[k, ]
kphat = phat[k, ]
if (kelimi[1] == 1 || sum(n[kelimi == 0]) ==
0) {
kseldose = 99
}
else {
adm.set = (kn != 0) & (kelimi == 0)
adm.index = which(adm.set == T)
y.adm = ky[adm.set]
n.adm = kn[adm.set]
selectd = sort(abs(kphat[adm.set] - target),
index.return = T)$ix[1]
kseldose = adm.index[selectd]
}
selectdoses[k, 2] = ifelse(is.na(kseldose),
99, kseldose)
if (k < dim(n)[1])
if (selectdoses[k + 1, 2] == dim(n)[2])
selectdoses[k, 2] = dim(n)[2]
if (k < dim(n)[1])
if (aa(selectdoses[k + 1, 2]) == dim(n)[2] &
aa(selectdoses[k + 1, 2]) == aa(selectdoses[k,
2]))
selectdoses[k, 2] = 99
}
}
else {
selectdoses = matrix(99, nrow = 1, ncol = 2)
selectdoses[1, ] = matrix(selectdose, nrow = 1)
}
selectdoses = matrix(selectdoses[selectdoses[, 2] !=
99, ], ncol = 2)
}
colnames(selectdoses) = c("DoseA", "DoseB")
}
if (mtd.contour == FALSE) {
if (selectdoses[1, 1] == 99 && selectdoses[1, 2] ==
99) {
out=list(target = target, MTD = 99, p_est = matrix(NA,
nrow = dim(npts)[1], ncol = dim(npts)[2]))
return(out)
}
else {
out=list(target = target, MTD = selectdoses,
p_est = round(phat.out, 2))
return(out)
}
}
else {
if (length(selectdoses) == 0) {
out=list(target = target, MTD = 99, p_est = matrix(NA,
nrow = dim(npts)[1], ncol = dim(npts)[2]))
return(out)
}
else {
out=list(target = target, MTD = selectdoses,
p_est = round(phat.out, 2))
return(out)
}
}
}
waterfall.subtrial.mtd <- function(target,npts, ntox, cutoff.eli = 0.95,
extrasafe = FALSE, offset = 0.05, boundMTD=FALSE, p.tox, temp) {
lambda_d = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -target)/(target * (1 - p.tox)))
b.e = temp[2, ]
pava <- function(x, wt = rep(1, length(x))) {
n <- length(x)
if (n <= 1)
return(x)
if (any(is.na(x)) || any(is.na(wt))) {
stop("Missing values in 'x' or 'wt' not allowed")
}
lvlsets <- (1:n)
repeat {
viol <- (as.vector(diff(x)) < 0)
if (!(any(viol)))
break
i <- min((1:(n - 1))[viol])
lvl1 <- lvlsets[i]
lvl2 <- lvlsets[i + 1]
ilvl <- (lvlsets == lvl1 | lvlsets == lvl2)
x[ilvl] <- sum(x[ilvl] * wt[ilvl])/sum(wt[ilvl])
lvlsets[ilvl] <- lvl1
}
x
}
y = ntox
n = npts
ndose = length(n)
elimi = rep(0, ndose)
is.escalation = 0
for (i in 1:ndose) {
if (n[i] >= 3) {
if (1 - pbeta(target, y[i] + 1, n[i] - y[i] +
1) > cutoff.eli) {
elimi[i:ndose] = 1
break
}
}
}
if (extrasafe) {
if (n[1] >= 3) {
if (1 - pbeta(target, y[1] + 1, n[1] - y[1] +
1) > cutoff.eli - offset) {
elimi[1:ndose] = 1
}
}
}
if (elimi[1] == 1 || sum(n[elimi == 0]) == 0) {
selectdose = 99
}
else {
adm.set = (n != 0) & (elimi == 0)
adm.index = which(adm.set == T)
y.adm = y[adm.set]
n.adm = n[adm.set]
phat = (y.adm + 0.05)/(n.adm + 0.1)
phat.var = (y.adm + 0.05) * (n.adm - y.adm + 0.05)/((n.adm +
0.1)^2 * (n.adm + 0.1 + 1))
phat = pava(phat, wt = 1/phat.var)
phat = phat + (1:length(phat)) * 1e-10
if(boundMTD){
if(all(phat>=lambda_d)){selectdose=99}else{
phat=phat[phat<lambda_d]
selectd = sort(abs(phat - target), index.return = T)$ix[1]
selectdose = adm.index[selectd]
}
}else{
selectd = sort(abs(phat - target), index.return = T)$ix[1]
selectdose = adm.index[selectd]
}
# selectd = sort(abs(phat - target), index.return = T)$ix[1]
# selectdose = adm.index[selectd]
if(selectdose!=99){
if (y[selectdose] <= b.e[n[selectdose]]) {
is.escalation = 1
}
}
}
list(selectdose = selectdose, is.escalation = is.escalation)
}
waterfall.subtrial <- function(target, p.true, dosespace,
npts, ntox, elimi, ncohort, cohortsize, n.earlystop = 20,
startdose = 1, p.saf = 0.6 * target, p.tox = 1.4 * target,
cutoff.eli = 0.95, extrasafe = FALSE, offset = 0.05,
totaln, titration.first.trial = FALSE, temp,boundMTD=FALSE) {
ndoses1 = nrow(p.true)
ndoses2 = ncol(p.true)
p.truee = p.true[dosespace]
npts = npts
ntox = ntox
elimi = elimi
if (target < 0.05) {
stop(" the target is too low! ")
}
if (target > 0.6) {
stop("the target is too high!")
}
if ((target - p.saf) < (0.1 * target)) {
stop("the probability deemed safe cannot be higher than or too close to the target! ")
}
if ((p.tox - target) < (0.1 * target)) {
stop("the probability deemed toxic cannot be lower than or too close to the target!")
}
if (offset >= 0.5) {
stop("the offset is too large! ")
}
ndose = length(p.truee)
selectdose = 0
is.escalation = 0
b.e = temp[2, ]
b.d = temp[3, ]
b.elim = temp[4, ]
lambda1 = log((1 - p.saf)/(1 - target))/log(target *
(1 - p.saf)/(p.saf * (1 - target)))
lambda2 = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -
target)/(target * (1 - p.tox)))
y <- rep(0, ndose)
n <- rep(0, ndose)
earlystop = 0
d = startdose
elm = rep(0, ndose)
if (titration.first.trial) {
z <- (runif(ndose) < p.truee)
if (sum(z) == 0) {
d = ndose
n[1:ndose] = 1
}
else {
d = which(z == 1)[1]
n[1:d] = 1
y[d] = 1
}
}
ft=TRUE #flag used to determine whether or not to add cohortsize-1 patients to a dose for the first time when titration is triggered.
for (icohort in 1:ncohort) {
if (titration.first.trial & n[d] < cohortsize & ft) {
ft=FALSE
y[d] = y[d] + sum(runif(cohortsize - 1) < p.truee[d])
n[d] = n[d] + cohortsize - 1
}else {
y[d] = y[d] + sum(runif(cohortsize) < p.truee[d])
n[d] = n[d] + cohortsize
}
if (!is.na(b.elim[n[d]])) {
if (y[d] >= b.elim[n[d]]) {
elm[d:ndose] = 1
if (d == 1) {
earlystop = 1
break
}
}
if (extrasafe) {
if (d == 1 && n[1] >= 3) {
if (1 - pbeta(target, y[1] + 0.05, n[1] -
y[1] + 0.1) > cutoff.eli - offset) {
earlystop = 1
break
}
}
}
}
if (y[d] <= b.e[n[d]] && d != ndose) {
if (elm[d + 1] == 0)
d = d + 1
}
else if (y[d] >= b.d[n[d]] && d != 1) {
d = d - 1
}
else {
d = d
}
if (n[d] >= n.earlystop)
break
if (sum(n) >= (ncohort * cohortsize))
break
}
if (earlystop == 1) {
selectdose = 99
elm = rep(1, ndose)
}
else {
wsmtd = waterfall.subtrial.mtd(target, n, y, cutoff.eli,
extrasafe, offset, boundMTD=boundMTD, p.tox,temp)
selectdose = wsmtd$selectdose
is.escalation = wsmtd$is.escalation
}
npts[dosespace] = n
ntox[dosespace] = y
elimi[dosespace] = elm
list(ncohort = icohort, ntotal = icohort * cohortsize,
startdose = startdose, npts = npts, ntox = ntox,
totaltox = sum(ntox), totaln = sum(npts), pctearlystop = sum(selectdose ==
99) * 100, selectdose = selectdose, is.escalation = is.escalation,
elimi = elimi)
}
get.oc.comb.waterfall <- function(p.true, target, ncohort,
cohortsize, n.earlystop = 12, cutoff.eli = 0.95, p.saf = 0.6 *
target, p.tox = 1.4 * target, titration = FALSE,
extrasafe = FALSE, offset = 0.05, boundMTD=FALSE,ntrial = 1000) {
temp = get.boundary(target, ncohort = 150, cohortsize = 1,
cutoff.eli = cutoff.eli, extrasafe = extrasafe)$boundary_tab
JJ = nrow(p.true)
KK = ncol(p.true)
if (JJ > KK)
p.true = t(p.true)
true.mtd.position = cbind(1:JJ, apply(p.true, 1, function(x) {
flaga = rep(0, length(x))
tmp = which.min(abs(x - target))
flaga[tmp] = 1
flagb = (x <= target + 0.05)
ifelse(sum(flaga & flagb) > 0, which(flaga & flagb),
99)
}))
true.mtd.pos = apply(true.mtd.position, 1, function(x) (x[2] -
1) * JJ + x[1])
true.mtd.pos.new = true.mtd.pos[true.mtd.pos <= JJ *
KK]
nMTDs = paste(sort(true.mtd.pos.new), collapse = ",")
greater.than.contour = rep(1, length(JJ * KK))
less.than.contour = 1 - greater.than.contour
if (sum(true.mtd.position[, 2] <= KK) > 0) {
tmp = true.mtd.position[true.mtd.position[, 2] <=
JJ * KK, ]
if (length(tmp) == 2)
tmp = matrix(tmp, ncol = 2)
ONES = matrix(1, nrow = JJ, ncol = KK)
for (kkk in 1:nrow(tmp)) {
tmpa = tmp[kkk, 1]
tmpb = tmp[kkk, 2]
ONES[1:tmpa, 1:tmpb] = 0
}
ONES[true.mtd.pos.new] = 1
less.than.contour = which(ONES == 0)
ONES = matrix(1, nrow = JJ, ncol = KK)
for (kkk in 1:nrow(tmp)) {
tmpa = tmp[kkk, 1]
tmpb = tmp[kkk, 2]
ONES[tmpa:JJ, tmpb:KK] = 0
}
ONES[true.mtd.pos.new] = 1
greater.than.contour = which(ONES == 0)
}
at.contour = c(1:(JJ * KK))[-c(greater.than.contour,
less.than.contour)]
aa = function(x) as.numeric(as.character(x))
if (target < 0.05) {
stop(" the target is too low! ")
}
if (target > 0.6) {
stop("the target is too high!")
}
if ((target - p.saf) < (0.1 * target)) {
stop("the probability deemed safe cannot be higher than or too close to the target! ")
}
if ((p.tox - target) < (0.1 * target)) {
stop("the probability deemed toxic cannot be lower than or too close to the target!")
}
ndoses1 <- nrow(p.true)
ndoses2 <- ncol(p.true)
ntrial.phase1 = NULL
ntrial.mtd = NULL
ntrial.nt = NULL
ntrial.yt = NULL
ntrial.ne = NULL
ntrial.ye = NULL
if (cohortsize == 1)
titration = FALSE
titration.first.trial = FALSE
for (trial in 1:ntrial) {
trial.result = NULL
ntox = matrix(rep(0, (ndoses2) * ndoses1), ncol = ndoses2)
colnames(ntox) = paste("ntoxDoseB", 1:ndoses2, sep = "")
npts = matrix(rep(0, (ndoses2) * ndoses1), ncol = ndoses2)
colnames(npts) = paste("nptsDoseB", 1:ndoses2, sep = "")
elimi = matrix(0, nrow = ndoses1, ncol = ndoses2)
colnames(elimi) = paste("elimiDoseB", 1:ndoses2,
sep = "")
mtd = cbind(selectdoseA = 1:ndoses1, selectdoseB = rep(NA,
ndoses1))
trial.result = data.frame(cbind(trial = rep(trial,
ndoses1), mtd, npts, ntox, elimi))
totaln = 0
startdose = 1
dosespace = c(1:(ndoses1 - 1), (1:ndoses2) * ndoses1)
subtriali = 1
while (totaln < sum(ncohort) * cohortsize) {
if (titration & subtriali == 1) {
titration.first.trial = TRUE
}
else {
titration.first.trial = FALSE
}
subtrial = waterfall.subtrial(target, p.true = p.true,
dosespace = dosespace, npts = npts, ntox = ntox,
elimi = elimi, ncohort = ncohort[subtriali],
cohortsize = cohortsize, n.earlystop = n.earlystop,
startdose = startdose, p.saf = p.saf, p.tox = p.tox,
cutoff.eli = cutoff.eli, extrasafe = extrasafe,
offset = offset, totaln = totaln,
titration.first.trial = titration.first.trial,
temp = temp,boundMTD=boundMTD)
selectdose = ifelse(subtrial$selectdose == 99,
99, dosespace[subtrial$selectdose])
if (selectdose == 99)
break
dj = ifelse(selectdose%%ndoses1 == 0, selectdose%/%ndoses1,
selectdose%/%ndoses1 + 1)
di = selectdose - (dj - 1) * ndoses1
totaln = aa(subtrial$totaln)
npts = subtrial$npts
ntox = subtrial$ntox
elimi = subtrial$elimi
if ((subtriali == 1) & (selectdose < ndoses1)) {
for (a in (di + 1):ndoses1) for (b in 1:ndoses2) elimi[a,
b] = 1
if (subtrial$is.escalation == 1) {
startdose = 1
dosespace1 = c(di + ((2:ndoses2) - 1) * ndoses1)
subtriali = subtriali + 1
subtrial1 = waterfall.subtrial(target, p.true = p.true,
dosespace = dosespace1, npts = npts, ntox = ntox,
elimi = elimi, ncohort = ncohort[subtriali],
cohortsize = cohortsize, n.earlystop = n.earlystop,
startdose = startdose, p.saf = p.saf, p.tox = p.tox,
cutoff.eli, extrasafe, offset, totaln = totaln,
temp = temp,boundMTD=boundMTD)
selectdose1 = ifelse(subtrial1$selectdose ==
99, selectdose, dosespace1[subtrial1$selectdose])
if (selectdose1 == 99)
break
dj = ifelse(selectdose1%%ndoses1 == 0, selectdose1%/%ndoses1,
selectdose1%/%ndoses1 + 1)
di = selectdose1 - (dj - 1) * ndoses1
totaln = aa(subtrial1$totaln)
npts = subtrial1$npts
ntox = subtrial1$ntox
elimi = subtrial1$elimi
}
}
if (di - 1 == 0)
break
subtriali = subtriali + 1
if (dj < ndoses2)
elimi[di, (dj + 1):ndoses2] = 1
startdose = dj
dosespace = di - 1 + ((2:ndoses2) - 1) * ndoses1
if (dj == ndoses2)
startdose = dj - 1
}
npts = t(apply(npts, 1, aa))
ntox = t(apply(ntox, 1, aa))
elimi = t(apply(elimi, 1, aa))
phat = (ntox + 0.05)/(npts + 0.1)
phat = t(apply(phat, 1, aa))
colnames(phat) = paste("phat", 1:ndoses2, sep = "")
phat[elimi == 1] = 1.1
phat = Iso::biviso(phat, npts + 0.1, warn = TRUE)[,
]
phat = phat + (1e-05) * (matrix(rep(1:dim(npts)[1],
each = dim(npts)[2], len = length(npts)), dim(npts)[1],
byrow = T) + matrix(rep(1:dim(npts)[2], each = dim(npts)[1],
len = length(npts)), dim(npts)[1]))
colnames(phat) = paste("phat", 1:ndoses2, sep = "")
for (k in ndoses1:1) {
kn = npts[k, ]
ky = ntox[k, ]
kelimi = elimi[k, ]
kphat = phat[k, ]
if (kelimi[1] == 1 || sum(npts[kelimi == 0]) ==
0) {
kseldose = 99
}
else {
adm.set = (kn != 0) & (kelimi == 0)
adm.index = which(adm.set == T)
y.adm = ky[adm.set]
n.adm = kn[adm.set]
selectd = sort(abs(kphat[adm.set] - target),
index.return = T)$ix[1]
kseldose = adm.index[selectd]
}
mtd[k, 2] = kseldose
if (k < ndoses1)
if (mtd[k, 2] - mtd[k + 1, 2] <= 0 & mtd[k +
1, 2] != 99)
mtd[k, 2] = mtd[k + 1, 2]
}
trial.result[1:ndoses1, grep("nptsDoseB", colnames(trial.result))] = npts
trial.result[1:ndoses1, grep("ntoxDoseB", colnames(trial.result))] = ntox
trial.result[1:ndoses1, grep("selectdose", colnames(trial.result))] = mtd
trial.result[1:ndoses1, grep("elimiDoseB", colnames(trial.result))] = elimi
ntrial.mtd = rbind(ntrial.mtd, cbind(trial = rep(trial,
nrow(mtd)), mtd))
ntrial.nt = rbind(ntrial.nt, cbind(trial = rep(trial,
nrow(npts)), npts))
ntrial.yt = rbind(ntrial.yt, cbind(trial = rep(trial,
nrow(ntox)), ntox))
trial.result = cbind(trial.result, phat)
ntrial.phase1 = rbind(ntrial.phase1, trial.result)
}
ntrial.mtd = data.frame(ntrial.mtd)
colnames(ntrial.mtd) = c("trial", "doseA", "doseB")
alltoxpercent = round(sum(sapply(1:ntrial, function(x) ifelse(sum(ntrial.mtd$doseB[ntrial.mtd$trial ==
x] == 99) == ndoses1, 1, 0)), na.rm = T) * 100/1000,
3)
stoppercent = round(sum(sapply(1:ntrial, function(x) sum(ntrial.mtd$doseB[ntrial.mtd$trial ==
x] == 99))) * 100/1000, 1)
selpercent = matrix(0, nrow = ndoses1, ncol = ndoses2)
nselpercent = 0
selpercent1 = 0
selpercent2 = 0
mtdtable = NULL
mtdlist = list()
for (triali in 1:ntrial) {
mtddata = unique(as.matrix(ntrial.mtd[ntrial.mtd$trial ==
triali, 2:3]))
mtddata = mtddata[!is.na(mtddata[, 2]) & mtddata[,
2] <= KK, ]
if (length(mtddata) > 0) {
if (length(mtddata) == 2)
mtddata = matrix(mtddata, ncol = 2)
mtdlevel = aa(t(apply(mtddata, 1, function(x) (x[2] -
1) * ndoses1 + x[1])))
mtdlevel[mtdlevel > ndoses1 * ndoses2] = 99
if (sum(mtdlevel <= ndoses1 * ndoses2) > 0) {
selpercent[mtdlevel[mtdlevel <= ndoses1 * ndoses2]] = selpercent[mtdlevel[mtdlevel <=
ndoses1 * ndoses2]] + 1
mtdlist[[triali]] = mtdlevel[mtdlevel <= ndoses1 *
ndoses2]
mtdtable = c(mtdtable, paste(sort(aa(mtdlevel[mtdlevel <=
ndoses1 * ndoses2])), collapse = ","))
if (paste(sort(aa(mtdlevel[mtdlevel <= ndoses1 *
ndoses2])), collapse = ",") == nMTDs)
nselpercent = nselpercent + 1
if (sum(mtdlevel <= ndoses1 * ndoses2) == 1 &
sum(is.element(mtdlevel[mtdlevel <= ndoses1 *
ndoses2], which(p.true == target)) == F) ==
0)
selpercent1 = selpercent1 + 1
if (sum(mtdlevel <= ndoses1 * ndoses2) == 2 &
sum(is.element(mtdlevel[mtdlevel <= ndoses1 *
ndoses2], which(p.true == target)) == F) ==
0)
selpercent2 = selpercent2 + 1
}
else {
mtdlist[[triali]] = 99
mtdtable = c(mtdtable, 99)
}
}
}
selpercent = round(selpercent * 100/ntrial, 2)
if (JJ > KK) {
nptsdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.nt), by = ndoses1)) nptsdose = nptsdose +
ntrial.nt[i + 0:(ndoses1 - 1), -1]
ntoxdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.yt), by = ndoses1)) ntoxdose = ntoxdose +
ntrial.yt[i + 0:(ndoses1 - 1), -1]
colnames(p.true) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(p.true) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(nptsdose) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(nptsdose) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(ntoxdose) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(ntoxdose) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
out = list(p.true = round(t(p.true), 2), selpercent = round(t(selpercent),
2), npatients = round(t(nptsdose)/ntrial, 2),
ntox = round(t(ntoxdose)/ntrial, 2), totaltox = round(sum(ntoxdose/ntrial),
1), totaln = round(sum(nptsdose/ntrial), 1))
out2 = list()
if (length(at.contour) > 0) {
out2 = list(npercent.contour = paste(round(100 *
sum(nptsdose[at.contour])/sum(nptsdose), 1),
"%", sep = ""), npercent.above.contour = paste(round(100 *
sum(nptsdose[greater.than.contour])/sum(nptsdose),
1), "%", sep = ""), npercent.below.contour = paste(round(100 *
sum(nptsdose[less.than.contour])/sum(nptsdose),
1), "%", sep = ""), pcs.contour = paste(round(100 *
nselpercent/ntrial, 1), "%", sep = ""), flowchart = FALSE)
}
outnew = c(out, out2)
class(outnew)<-"boin"
return(outnew)
}
else {
nptsdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.nt), by = ndoses1)) nptsdose = nptsdose +
ntrial.nt[i + 0:(ndoses1 - 1), -1]
ntoxdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.yt), by = ndoses1)) ntoxdose = ntoxdose +
ntrial.yt[i + 0:(ndoses1 - 1), -1]
out2 = list()
if (length(at.contour) > 0) {
out2 = list(npercent.contour = paste(round(100 *
sum(nptsdose[at.contour])/sum(nptsdose), 1),
"%", sep = ""), npercent.above.contour = paste(round(100 *
sum(nptsdose[greater.than.contour])/sum(nptsdose),
1), "%", sep = ""), npercent.below.contour = paste(round(100 *
sum(nptsdose[less.than.contour])/sum(nptsdose),
1), "%", sep = ""), pcs.contour = paste(round(100 *
nselpercent/ntrial, 1), "%", sep = ""))
}
rownames(p.true) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(p.true) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(nptsdose) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(nptsdose) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(ntoxdose) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(ntoxdose) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
out = list(p.true = apply(formatC(p.true, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
selpercent = apply(formatC(selpercent, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
npatients = apply(formatC(nptsdose/ntrial, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
ntox = apply(formatC(ntoxdose/ntrial, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
totaltox = as.numeric(formatC(sum(ntoxdose/ntrial),
digits = 1, format = "f")), totaln = as.numeric(formatC(sum(nptsdose/ntrial),
digits = 1, format = "f")), flowchart = FALSE)
outnew = c(out, out2)
class(outnew)<-"boin"
return(outnew)
}
}
set.seed(seed)
JJ = nrow(p.true)
KK = ncol(p.true)
if (JJ > KK) {
stop("p.true should be arranged in a way (i.e., rotated) such that the number of rows is less than or equal to the number of columns.")
}
if (mtd.contour == FALSE) {
if (is.null(n.earlystop) == TRUE)
n.earlystop = 100
if (n.earlystop <= 6) {
warning("the value of n.earlystop is too low to ensure good operating characteristics. Recommend n.earlystop = 9 to 18 ")
}
if (length(ncohort) > 1) {
warning("ncohort is the total number of cohorts for the trial. Please enter a scalar.")
}
if (((JJ * KK) <= 4) & (sum(ncohort) <= 6)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
if (((JJ * KK) > 4) & (sum(ncohort) <= 8)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
out=get.oc.comb.boin(target = target, p.true = p.true,
ncohort = sum(ncohort), cohortsize = cohortsize,
n.earlystop = n.earlystop, startdose = startdose,
titration = titration,
p.saf = p.saf, p.tox = p.tox,
cutoff.eli = cutoff.eli, extrasafe = extrasafe,
offset = offset,boundMTD=boundMTD,
ntrial = ntrial)
class(out)<-"boin"
return(out)
}
if (mtd.contour == TRUE) {
if (missing(ncohort) == TRUE) {
constSeq = rep(round(4/cohortsize, 2), 20)
dosespaceSeq = c(JJ + KK - 1, rep(KK - 1, JJ - 1))
ncohort = ceiling(constSeq[1:JJ] * dosespaceSeq)
}
if (length(ncohort) != JJ) {
stop("The vector length of ncohort doesn't match the number of subtrials (the number of dose matrix rows).\n Please enter the number of ncohorts of each subtrial.")
}
if (((JJ * KK) <= 4) & (sum(ncohort) <= 6)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
if (((JJ * KK) > 4) & (sum(ncohort) <= 8)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
if (is.null(n.earlystop) == TRUE)
n.earlystop = 12
if (n.earlystop <= 6) {
warning("the value of n.earlystop is too low to ensure good operating characteristics. Recommend n.earlystop = 9 to 18 ")
}
out=get.oc.comb.waterfall(p.true = p.true, target = target,
ncohort = ncohort, cohortsize = cohortsize, n.earlystop = n.earlystop,
cutoff.eli = cutoff.eli, p.saf = p.saf, p.tox = p.tox,
titration = titration,
extrasafe = extrasafe, offset = offset,boundMTD=boundMTD,
ntrial = ntrial)
class(out)<-"boin"
return(out)
}
}
Try the BOIN package in your browser
Any scripts or data that you put into this service are public.
BOIN documentation built on Jan. 20, 2021, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#'
#' Generate operating characteristics for drug combination trials
#'
#' Obtain the operating characteristics of the BOIN design or waterfall design for drug combination
#' trials. The BOIN design is to find a MTD, and the waterfall design is to find the MTD contour
#' (i.e., multple MTDs in the dose matrix)
#'
#' @usage get.oc.comb(target, p.true, ncohort, cohortsize, n.earlystop=NULL, startdose=c(1, 1),
#' titration=FALSE,p.saf=0.6*target, p.tox=1.4*target, cutoff.eli=0.95,
#' extrasafe=FALSE,offset=0.05, ntrial=1000, mtd.contour=FALSE,
#' boundMTD=FALSE, seed=6)
#'
#' @param target the target DLT rate
#' @param p.true a \code{J*K} matrix \code{(J<=K)} containing the true toxicity probabilities of
#' combinations with \code{J} dose levels of agent A and \code{K} dose levels of agent B
#' @param ncohort a \code{1*J} vector specifying the number of cohorts for each of \code{J} subtrials
#' if \code{mtd.contour=TRUE}; Otherwise, a scalar specifying the total number of cohorts for
#' the trial.
#' @param cohortsize the cohort size
#' @param n.earlystop the early stopping parameter. If the number of patients treated at the current
#' dose reaches \code{n.earlystop}, stop the trial or subtrial and select the MTD based on
#' the observed data. When the waterfall design is used to find the MTD contour,
#' \code{n.earlystop=12} by default.
#' @param startdose the starting dose combination level for drug combination trial
#' @param titration set \code{titration=TRUE} to perform dose escalation with cohort size = 1 to accelerate dose escalation at the begining of the trial.
#' @param p.saf the highest toxicity probability that is deemed subtherapeutic (i.e. below the MTD)
#' such that dose escalation should be undertaken.
#' The default value is \code{p.saf=0.6*target}.
#' @param p.tox the lowest toxicity probability that is deemed overly toxic such that deescalation
#' is required. The default value is \code{p.tox=1.4*target}.
#' @param cutoff.eli the cutoff to eliminate an overly toxic dose for safety. We recommend the
#' default value of (\code{cutoff.eli=0.95}) for general use.
#' @param extrasafe set \code{extrasafe=TRUE} to impose a more stringent stopping rule
#' @param offset a small positive number (between 0 and 0.5) to control how strict the stopping
#' rule is when \code{extrasafe=TRUE}. A larger value leads to a more strict stopping
#' rule. The default value \code{offset=0.05} generally works well.
#' @param ntrial the total number of trials to be simulated
#' @param mtd.contour set \code{mtd.contour=TRUE} to select the MTD contour (claiming multiple MTDs).
#' Otherwise, BOIN design is used to search for a single MTD.
#' @param boundMTD set \code{boundMTD=TRUE} to impose the condition: the isotonic estimate of toxicity
#' probability for the selected MTD must be less than de-escalation boundary.
#' @param seed the random seed for simulation
#'
#' @details The operating characteristics of the BOIN design or waterfall design are generated by
#' simulating trials under the prespecified true toxicity probabilities of the investigational dose
#' combinations. If \code{titration=TRUE}, we perform dose escalation with cohort size = 1 at the begining of the trial:
#' starting from \code{startdose}, if no toxicity is observed, we escalate the dose;
#' otherwise, the titration is completed and we switch to cohort size = \code{cohortsize}.
#' Titration accelerates the dose escalation and is useful when low doses are believed to be safe.
#'
#'
#' The BOIN and waterfall designs have two built-in stopping rules:
#' (1) stop the trial/subtrial if the lowest dose is eliminated due to toxicity, and no dose should
#' be selected as the MTD; and (2) stop the trial/subtrial and select the MTD if the number of
#' patients treated at the current dose reaches \code{n.earlystop}. The first stopping rule is a safety
#' rule to protect patients from the case in which all doses are overly toxic. The rationale for
#' the second stopping rule is that when there is a large number (i.e., \code{n.earlystop}) of
#' patients assigned to a dose, it means that the dose-finding algorithm has approximately converged.
#' Thus, we can stop the trial/subtrial early and select the MTD to save sample size and reduce the
#' trial duration.
#'
#'
#' For some applications, investigators may prefer a more strict safety stopping rule than rule
#' (1) for extra safety when the lowest dose is overly toxic.
#' This can be achieved by setting \code{extrasafe=TRUE},
#' which imposes the following more strict safety stopping rule:
#' stop the trial if (i) the number of patients treated at the lowest dose \code{>=3},
#' and (ii) \eqn{Pr(toxicity\ rate\ of\ the\ lowest\ dose > \code{target} | data) > \code{cutoff.eli}-\code{offset}}.
#' As a tradeoff, the strong stopping rule will decrease the MTD selection percentage
#' when the lowest dose actually is the MTD.
#'
#' @return \code{get.oc.comb()} returns the operating characteristics of the BOIN combination or
#' waterfall design as a list. For the BOIN combination design, including:
#' (1) true toxicity probability at each dose level (\code{$p.true}),
#' (2) selection percentage at each dose level (\code{$selpercent}),
#' (3) the number of patients treated at each dose level (\code{$npatients})
#' (4) the number of toxicities observed at each dose level (\code{$ntox}),
#' (5) the total number of patients in the trial (\code{$totaln}),
#' (6) the total number of toxicities observed for the trial (\code{$totaltox})
#' (7) the pecentage of correct selection (\code{$pcs}),
#' (8) the total percentage of patients treated at the MTD (\code{$npercent}).
#' (9) the percentage of early stopping without selecting the MTD (\code{$percentstop})
#' For the the waterfall design, including:
#' (1) true toxicity probability at each dose level (\code{$p.true}),
#' (2) selection percentage of dose combinations (\code{$selpercent}),
#' (3) the number of patients treated at each dose combination (\code{$npatients})
#' (4) the number of toxicities observed at each dose combination (\code{$ntox}),
#' (5) the total number of patients in the trial (\code{$totaln}),
#' (6) the total number of toxicities observed for the trial (\code{$totaltox})
#' (7) the total percentage of correct selection at the MTD contour (\code{$pcs.contour}),
#' (8) the total percentage of patients treated at MTD contour
#' (\code{$npercent.contour})
#' (9) the total percentage of patients treated above MTD contour
#' (\code{$npercent.above.contour})
#' (10) the total percentage of patients treated below MTD contour
#' (\code{$npercent.below.contour})
#'
#'
#' @note We should avoid setting the values of \code{p.saf} and \code{p.tox} very close to the
#' \code{target}. This is because the small sample sizes of typical phase I trials prevent us from
#' differentiating the target DLT rate from the rates close to it. The default values provided by
#' \code{get.oc()} are strongly recommended, and generally yield excellent operating characteristics.
#'
#'
#' @author Suyu Liu, Liangcai Zhang, Yanhong Zhou, and Ying Yuan
#'
#' @references Liu S. and Yuan, Y. (2015). Bayesian Optimal Interval Designs for Phase I Clinical
#' Trials, \emph{Journal of the Royal Statistical Society: Series C}, 64, 507-523.
#'
#' Lin R. and Yin, G. (2017). Bayesian Optimal Interval Designs for Dose Finding in
#' Drug-combination Trials, \emph{Statistical Methods in Medical Research}, 26, 2155-2167.
#'
#' Yan, F., Zhang, L., Zhou, Y., Pan, H., Liu, S. and Yuan, Y. (2020).BOIN: An R Package
#' for Designing Single-Agent and Drug-Combination Dose-Finding Trials Using Bayesian Optimal
#' Interval Designs. \emph{Journal of Statistical Software}, 94(13),1-32.<doi:10.18637/jss.v094.i13>.
#'
#'
#' Zhang L. and Yuan, Y. (2016). A Simple Bayesian Design to Identify the Maximum
#' Tolerated Dose Contour for Drug Combination Trials, \emph{Statistics in Medicine}, 35, 4924-4936.
#'
#' @seealso Tutorial: \url{http://odin.mdacc.tmc.edu/~yyuan/Software/BOIN/BOIN2.6_tutorial.pdf}
#'
#' Paper: \url{http://odin.mdacc.tmc.edu/~yyuan/Software/BOIN/paper.pdf}
#'
#' @examples
#'
#' ###### drug-combination trial ######
#'
#' ##### combination trial to find a single MTD ######
#'
#' ## get the operating characteristics for BOIN design
#' p.true <- matrix(c(0.01,0.03,0.10,0.20,0.30,
#' 0.03,0.05,0.15,0.30,0.60,
#' 0.08,0.10,0.30,0.60,0.75), byrow=TRUE, ncol=5)
#'
#' oc.comb <- get.oc.comb(target=0.3, p.true, ncohort=20, cohortsize=3,
#' n.earlystop=12, startdose=c(1,1), ntrial=100)
#'
#' summary(oc.comb)
#' plot(oc.comb)
#'
#'
#' ## get the operating characteristics with titration for BOIN design
#' oc.comb <- get.oc.comb(target=0.3, p.true, ncohort=20, cohortsize=3,
#' n.earlystop=12, startdose=c(1,1), titration=TRUE, ntrial=100)
#' summary(oc.comb)
#' plot(oc.comb)
#'
#'
#' ##### combination trial to find the MTD contour ######
#'
#' ## find the MTD contour using waterfall design
#' oc.comb <- get.oc.comb(target=0.3, p.true, ncohort=c(10,5,5), cohortsize=3,
#' n.earlystop=12, startdose=c(1,1), ntrial=100, mtd.contour=TRUE)
#'
#' summary(oc.comb)
#' plot(oc.comb)
#'
#' @export
get.oc.comb <- function (target, p.true, ncohort, cohortsize, n.earlystop = NULL,
startdose = c(1, 1), titration = FALSE, p.saf = 0.6 * target,
p.tox = 1.4 * target, cutoff.eli = 0.95, extrasafe = FALSE,
offset = 0.05, ntrial = 1000, mtd.contour = FALSE, boundMTD=FALSE,seed = 6)
{
get.oc.comb.boin <- function(target, p.true, ncohort, cohortsize,
n.earlystop = 100, startdose = c(1, 1), titration = FALSE,
p.saf = 0.6 * target, p.tox = 1.4 * target, cutoff.eli = 0.95,
extrasafe = FALSE, offset = 0.05, boundMTD=FALSE,ntrial = 1000) {
JJ = nrow(p.true)
KK = ncol(p.true)
if (JJ > KK)
stop("p.true should be arranged in a way (i.e., rotated) such that\n the number of rows is less than or equal to the number of columns.")
if (JJ > KK)
p.true = t(p.true)
if (target < 0.05) {
stop("the target is too low!")
}
if (target > 0.6) {
stop("the target is too high!")
}
if ((target - p.saf) < (0.1 * target)) {
stop("the probability deemed safe cannot be higher than or too close to the target!")
}
if ((p.tox - target) < (0.1 * target)) {
stop("the probability deemed toxic cannot be lower than or too close to the target!")
}
if (offset >= 0.5) {
stop("the offset is too large!")
}
if (n.earlystop <= 6) {
warning("the value of n.earlystop is too low to ensure good operating characteristics. Recommend n.earlystop = 9 to 18")
}
ndose = length(p.true)
npts = ncohort * cohortsize
Y <- array(matrix(rep(0, length(p.true) * ntrial), dim(p.true)[1]),
dim = c(dim(p.true), ntrial))
N <- array(matrix(rep(0, length(p.true) * ntrial), dim(p.true)[1]),
dim = c(dim(p.true), ntrial))
dselect = matrix(rep(0, 2 * ntrial), ncol = 2)
if (cohortsize > 1) {
temp = get.boundary(target, ncohort, cohortsize,
n.earlystop=100, p.saf, p.tox, cutoff.eli, extrasafe)$full_boundary_tab
}else {
temp = get.boundary(target, ncohort, cohortsize,
n.earlystop=100, p.saf, p.tox, cutoff.eli, extrasafe)$boundary_tab
}
b.e = temp[2, ]
b.d = temp[3, ]
b.elim = temp[4, ]
lambda1 = log((1 - p.saf)/(1 - target))/log(target *
(1 - p.saf)/(p.saf * (1 - target)))
lambda2 = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -
target)/(target * (1 - p.tox)))
if (cohortsize == 1)
titration = FALSE
for (trial in 1:ntrial) {
y <- matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
n <- matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
earlystop = 0
d = startdose
elimi = matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
ft=TRUE #flag used to determine whether or not to add cohortsize-1 patients to a dose for the first time when titration is triggered.
if (titration) {
tmpa = d[1]
tmpb = d[2]
y[tmpa, tmpb] <- (runif(1) < p.true[tmpa, tmpb])
n[tmpa, tmpb] <- 1
while (tmpa <= dim(p.true)[1] & tmpb <= dim(p.true)[2]) {
if (tmpa == dim(p.true)[1] & tmpb == dim(p.true)[2]) {
break
}
if (sum(y) == 1) {
y[tmpa, tmpb] = 1
break
}
if (tmpa < dim(p.true)[1] & tmpb < dim(p.true)[2]) {
tmp.candidate = rbind(c(tmpa + 1, tmpb),
c(tmpa, tmpb + 1))
tmp.sel = rbinom(1, 1, prob = c(0.5, 0.5)) +
1
tmpa = tmp.candidate[tmp.sel, 1]
tmpb = tmp.candidate[tmp.sel, 2]
}
else if (tmpa == dim(p.true)[1]) {
tmpb = tmpb + 1
}
else {
tmpa = tmpa + 1
}
y[tmpa, tmpb] <- (runif(1) < p.true[tmpa, tmpb])
n[tmpa, tmpb] <- 1
}
if (sum(y) == 0) {
d = c(dim(p.true)[1], dim(p.true)[2])
}
else {
d = c(tmpa, tmpb)
}
}
for (pp in 1:ncohort) {
if (titration & n[d[1], d[2]] < cohortsize & ft) {
ft=FALSE
y[d[1], d[2]] = y[d[1], d[2]] + sum(runif(cohortsize -
1) < p.true[d[1], d[2]])
n[d[1], d[2]] = n[d[1], d[2]] + cohortsize -1
}
else {
y[d[1], d[2]] = y[d[1], d[2]] + sum(runif(cohortsize) <
p.true[d[1], d[2]])
n[d[1], d[2]] = n[d[1], d[2]] + cohortsize
}
nc = n[d[1], d[2]]
if (!is.na(b.elim[nc])) {
if (y[d[1], d[2]] >= b.elim[nc]) {
for (i in min(d[1], dim(p.true)[1]):dim(p.true)[1]) {
for (j in min(d[2], dim(p.true)[2]):dim(p.true)[2]) {
elimi[i, j] = 1
}
}
if (d[1] == 1 && d[2] == 1) {
d = c(99, 99)
earlystop = 1
break
}
}
if (extrasafe) {
if (d[1] == 1 && d[2] == 1 && n[1, 1] >=
3) {
if (1 - pbeta(target, y[1, 1] + 1, n[1,
1] - y[1, 1] + 1) > cutoff.eli - offset) {
d = c(99, 99)
earlystop = 1
break
}
}
}
}
if (n[d[1],d[2]] >= n.earlystop && (y[d[1],d[2]]>b.e[n[d[1],d[2]]] ||
(d[1]==dim(p.true)[1] && d[2]==dim(p.true)[2]) ||
( d[1]==dim(p.true)[1] && d[2]<dim(p.true)[2] && elimi[d[1],d[2]+1]==1 ) ||
( d[1]<dim(p.true)[1] && d[2]==dim(p.true)[2] && elimi[d[1]+1,d[2]]==1 ) ||
( d[1]<dim(p.true)[1] && d[2]<dim(p.true)[2] && elimi[d[1]+1,d[2]]==1 && elimi[d[1],d[2]+1]==1 ) ) &&
(y[d[1],d[2]]<b.d[n[d[1],d[2]]] || (d[1]==1 && d[2]==1) ) ) break;
if (y[d[1], d[2]] <= b.e[nc]) {
elevel = matrix(c(1, 0, 0, 1), 2)
pr_H0 = rep(0, length(elevel)/2)
nn = pr_H0
for (i in seq(1, length(elevel)/2, by = 1)) {
if (d[1] + elevel[1, i] <= dim(p.true)[1] &&
d[2] + elevel[2, i] <= dim(p.true)[2]) {
if (elimi[d[1] + elevel[1, i], d[2] + elevel[2,
i]] == 0) {
yn = y[d[1] + elevel[1, i], d[2] + elevel[2,
i]]
nn[i] = n[d[1] + elevel[1, i], d[2] +
elevel[2, i]]
pr_H0[i] <- pbeta(lambda2, yn + 0.5,
nn[i] - yn + 0.5) - pbeta(lambda1,
yn + 0.5, nn[i] - yn + 0.5)
}
}
}
pr_H0 = pr_H0 + nn * 5e-04
if (max(pr_H0) == 0) {
d = d
}
else {
k = which(pr_H0 == max(pr_H0))[as.integer(runif(1) *
length(which(pr_H0 == max(pr_H0))) + 1)]
d = d + c(elevel[1, k], elevel[2, k])
}
}
else if (y[d[1], d[2]] >= b.d[nc]) {
delevel = matrix(c(-1, 0, 0, -1), 2)
pr_H0 = rep(0, length(delevel)/2)
nn = pr_H0
for (i in seq(1, length(delevel)/2, by = 1)) {
if (d[1] + delevel[1, i] > 0 && d[2] + delevel[2,
i] > 0) {
yn = y[d[1] + delevel[1, i], d[2] + delevel[2,
i]]
nn[i] = n[d[1] + delevel[1, i], d[2] +
delevel[2, i]]
pr_H0[i] = pbeta(lambda2, yn + 0.5, nn[i] -
yn + 0.5) - pbeta(lambda1, yn + 0.5,
nn[i] - yn + 0.5)
}
}
pr_H0 = pr_H0 + nn * 5e-04
if (max(pr_H0) == 0) {
d = d
}
else {
k = which(pr_H0 == max(pr_H0))[as.integer(runif(1) *
length(which(pr_H0 == max(pr_H0))) + 1)]
d = d + c(delevel[1, k], delevel[2, k])
}
}
else {
d = d
}
}
Y[, , trial] = y
N[, , trial] = n
if (earlystop == 1) {
dselect[trial, ] = c(99, 99)
}else {
selcomb = select.mtd.comb.boin(target, n, y,
cutoff.eli, extrasafe, offset,
boundMTD=boundMTD,p.tox=p.tox,mtd.contour = FALSE)$MTD
dselect[trial, 1] = selcomb[1]
dselect[trial, 2] = selcomb[2]
}
}
selpercent = matrix(rep(0, ndose), dim(p.true)[1], dim(p.true)[2])
nptsdose = apply(N, c(1, 2), mean, digits = 2, format = "f")
ntoxdose = apply(Y, c(1, 2), mean, digits = 2, format = "f")
for (i in 1:dim(p.true)[1]) for (j in 1:dim(p.true)[2]) {
{
selpercent[i, j] = sum(dselect[, 1] == i & dselect[,
2] == j)/ntrial * 100
}
}
if (JJ <= KK) {
rownames(p.true) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(p.true) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
out = list(p.true = round(p.true, 2), selpercent = round(selpercent,2),
npatients = round(apply(N, c(1, 2), mean),2), ntox = round(apply(Y, c(1, 2), mean), 2),
totaltox = round(sum(Y)/ntrial, 1), totaln = round(sum(N)/ntrial,1),
pcs = paste(round(sum(selpercent[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)]),1), "%", sep = ""),
npercent = paste(round(sum(nptsdose[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)])/sum(nptsdose) *100, 1), "%", sep = ""),
percentstop=100-sum(round(selpercent,2)),flowchart = FALSE)
rownames(out$npatients) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(out$npatients) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(out$ntox) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(out$ntox) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
return(out)
}
else {
colnames(p.true) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(p.true) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(npatients) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(npatients) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(ntox) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(ntox) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
out = list(p.true = round(t(p.true), 2), selpercent = round(t(selpercent),2),
npatients = round(t(apply(N, c(1, 2), mean)), 2), ntox = round(t(apply(Y, c(1, 2), mean)),2),
totaltox = round(sum(Y)/ntrial, 1), totaln = round(sum(N)/ntrial,1), pcs = paste(round(sum(selpercent[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)]),
1), "%"), npercent = paste(round(sum(nptsdose[which(abs(p.true -target) == min(abs(p.true - target)), arr.ind = TRUE)])/sum(nptsdose) *100, 1), "%"),
percentstop=100-sum(round(selpercent,2)),
flowchart = FALSE)
return(out)
}
}
select.mtd.comb.boin <- function(target, npts, ntox, cutoff.eli = 0.95,
extrasafe = FALSE, offset = 0.05,
boundMTD=FALSE,p.tox=1.4*target,
mtd.contour = FALSE) {
lambda_d = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -target)/(target * (1 - p.tox)))
y = ntox
n = npts
if (nrow(n) > ncol(n) | nrow(y) > ncol(y)) {
stop("npts and ntox should be arranged in a way (i.e., rotated) such that for each of them, the number of rows is less than or equal to the number of columns.")
}
elimi = matrix(0, dim(n)[1], dim(n)[2])
if (extrasafe) {
if (n[1, 1] >= 3) {
if (1 - pbeta(target, y[1, 1] + 1, n[1, 1] -
y[1, 1] + 1) > cutoff.eli - offset) {
elimi[, ] = 1
}
}
}
for (i in 1:dim(n)[1]) {
for (j in 1:dim(n)[2]) {
if (n[i, j] >= 3) {
if (1 - pbeta(target, y[i, j] + 1, n[i, j] -
y[i, j] + 1) > cutoff.eli) {
elimi[i:dim(n)[1], j] = 1
elimi[i, j:dim(n)[2]] = 1
break
}
}
}
}
selectdose=NULL
if (elimi[1] == 1) {
selectdose = c(99, 99)
selectdoses = matrix(selectdose, nrow = 1)
}else {
phat = (y + 0.05)/(n + 0.1)
phat = Iso::biviso(phat, n + 0.1, warn = TRUE)[,
]
phat.out = phat
phat.out[n == 0] = NA
phat[elimi == 1] = 1.1
phat = phat * (n != 0) + (1e-05) * (matrix(rep(1:dim(n)[1],
each = dim(n)[2], len = length(n)), dim(n)[1],
byrow = T) + matrix(rep(1:dim(n)[2], each = dim(n)[1],
len = length(n)), dim(n)[1]))
if(boundMTD){
if(all(phat[n!=0]>=lambda_d)){
selectdose = c(99, 99)
selectdoses = matrix(selectdose, nrow = 1)
}else{
phat[phat>=lambda_d]=10}}
if(is.null(selectdose)){
phat[n == 0] = 10
selectdose = which(abs(phat - target) == min(abs(phat -
target)), arr.ind = TRUE)
if (length(selectdose) > 2)
selectdose = selectdose[1, ]
aa = function(x) as.numeric(as.character(x))
if (mtd.contour == TRUE) {
selectdoses = cbind(row = 1:dim(n)[1], col = rep(99,
dim(n)[1]))
for (k in dim(n)[1]:1) {
kn = n[k, ]
ky = y[k, ]
kelimi = elimi[k, ]
kphat = phat[k, ]
if (kelimi[1] == 1 || sum(n[kelimi == 0]) ==
0) {
kseldose = 99
}
else {
adm.set = (kn != 0) & (kelimi == 0)
adm.index = which(adm.set == T)
y.adm = ky[adm.set]
n.adm = kn[adm.set]
selectd = sort(abs(kphat[adm.set] - target),
index.return = T)$ix[1]
kseldose = adm.index[selectd]
}
selectdoses[k, 2] = ifelse(is.na(kseldose),
99, kseldose)
if (k < dim(n)[1])
if (selectdoses[k + 1, 2] == dim(n)[2])
selectdoses[k, 2] = dim(n)[2]
if (k < dim(n)[1])
if (aa(selectdoses[k + 1, 2]) == dim(n)[2] &
aa(selectdoses[k + 1, 2]) == aa(selectdoses[k,
2]))
selectdoses[k, 2] = 99
}
}
else {
selectdoses = matrix(99, nrow = 1, ncol = 2)
selectdoses[1, ] = matrix(selectdose, nrow = 1)
}
selectdoses = matrix(selectdoses[selectdoses[, 2] !=
99, ], ncol = 2)
}
colnames(selectdoses) = c("DoseA", "DoseB")
}
if (mtd.contour == FALSE) {
if (selectdoses[1, 1] == 99 && selectdoses[1, 2] ==
99) {
out=list(target = target, MTD = 99, p_est = matrix(NA,
nrow = dim(npts)[1], ncol = dim(npts)[2]))
return(out)
}
else {
out=list(target = target, MTD = selectdoses,
p_est = round(phat.out, 2))
return(out)
}
}
else {
if (length(selectdoses) == 0) {
out=list(target = target, MTD = 99, p_est = matrix(NA,
nrow = dim(npts)[1], ncol = dim(npts)[2]))
return(out)
}
else {
out=list(target = target, MTD = selectdoses,
p_est = round(phat.out, 2))
return(out)
}
}
}
waterfall.subtrial.mtd <- function(target,npts, ntox, cutoff.eli = 0.95,
extrasafe = FALSE, offset = 0.05, boundMTD=FALSE, p.tox, temp) {
lambda_d = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -target)/(target * (1 - p.tox)))
b.e = temp[2, ]
pava <- function(x, wt = rep(1, length(x))) {
n <- length(x)
if (n <= 1)
return(x)
if (any(is.na(x)) || any(is.na(wt))) {
stop("Missing values in 'x' or 'wt' not allowed")
}
lvlsets <- (1:n)
repeat {
viol <- (as.vector(diff(x)) < 0)
if (!(any(viol)))
break
i <- min((1:(n - 1))[viol])
lvl1 <- lvlsets[i]
lvl2 <- lvlsets[i + 1]
ilvl <- (lvlsets == lvl1 | lvlsets == lvl2)
x[ilvl] <- sum(x[ilvl] * wt[ilvl])/sum(wt[ilvl])
lvlsets[ilvl] <- lvl1
}
x
}
y = ntox
n = npts
ndose = length(n)
elimi = rep(0, ndose)
is.escalation = 0
for (i in 1:ndose) {
if (n[i] >= 3) {
if (1 - pbeta(target, y[i] + 1, n[i] - y[i] +
1) > cutoff.eli) {
elimi[i:ndose] = 1
break
}
}
}
if (extrasafe) {
if (n[1] >= 3) {
if (1 - pbeta(target, y[1] + 1, n[1] - y[1] +
1) > cutoff.eli - offset) {
elimi[1:ndose] = 1
}
}
}
if (elimi[1] == 1 || sum(n[elimi == 0]) == 0) {
selectdose = 99
}
else {
adm.set = (n != 0) & (elimi == 0)
adm.index = which(adm.set == T)
y.adm = y[adm.set]
n.adm = n[adm.set]
phat = (y.adm + 0.05)/(n.adm + 0.1)
phat.var = (y.adm + 0.05) * (n.adm - y.adm + 0.05)/((n.adm +
0.1)^2 * (n.adm + 0.1 + 1))
phat = pava(phat, wt = 1/phat.var)
phat = phat + (1:length(phat)) * 1e-10
if(boundMTD){
if(all(phat>=lambda_d)){selectdose=99}else{
phat=phat[phat<lambda_d]
selectd = sort(abs(phat - target), index.return = T)$ix[1]
selectdose = adm.index[selectd]
}
}else{
selectd = sort(abs(phat - target), index.return = T)$ix[1]
selectdose = adm.index[selectd]
}
# selectd = sort(abs(phat - target), index.return = T)$ix[1]
# selectdose = adm.index[selectd]
if(selectdose!=99){
if (y[selectdose] <= b.e[n[selectdose]]) {
is.escalation = 1
}
}
}
list(selectdose = selectdose, is.escalation = is.escalation)
}
waterfall.subtrial <- function(target, p.true, dosespace,
npts, ntox, elimi, ncohort, cohortsize, n.earlystop = 20,
startdose = 1, p.saf = 0.6 * target, p.tox = 1.4 * target,
cutoff.eli = 0.95, extrasafe = FALSE, offset = 0.05,
totaln, titration.first.trial = FALSE, temp,boundMTD=FALSE) {
ndoses1 = nrow(p.true)
ndoses2 = ncol(p.true)
p.truee = p.true[dosespace]
npts = npts
ntox = ntox
elimi = elimi
if (target < 0.05) {
stop(" the target is too low! ")
}
if (target > 0.6) {
stop("the target is too high!")
}
if ((target - p.saf) < (0.1 * target)) {
stop("the probability deemed safe cannot be higher than or too close to the target! ")
}
if ((p.tox - target) < (0.1 * target)) {
stop("the probability deemed toxic cannot be lower than or too close to the target!")
}
if (offset >= 0.5) {
stop("the offset is too large! ")
}
ndose = length(p.truee)
selectdose = 0
is.escalation = 0
b.e = temp[2, ]
b.d = temp[3, ]
b.elim = temp[4, ]
lambda1 = log((1 - p.saf)/(1 - target))/log(target *
(1 - p.saf)/(p.saf * (1 - target)))
lambda2 = log((1 - target)/(1 - p.tox))/log(p.tox * (1 -
target)/(target * (1 - p.tox)))
y <- rep(0, ndose)
n <- rep(0, ndose)
earlystop = 0
d = startdose
elm = rep(0, ndose)
if (titration.first.trial) {
z <- (runif(ndose) < p.truee)
if (sum(z) == 0) {
d = ndose
n[1:ndose] = 1
}
else {
d = which(z == 1)[1]
n[1:d] = 1
y[d] = 1
}
}
ft=TRUE #flag used to determine whether or not to add cohortsize-1 patients to a dose for the first time when titration is triggered.
for (icohort in 1:ncohort) {
if (titration.first.trial & n[d] < cohortsize & ft) {
ft=FALSE
y[d] = y[d] + sum(runif(cohortsize - 1) < p.truee[d])
n[d] = n[d] + cohortsize - 1
}else {
y[d] = y[d] + sum(runif(cohortsize) < p.truee[d])
n[d] = n[d] + cohortsize
}
if (!is.na(b.elim[n[d]])) {
if (y[d] >= b.elim[n[d]]) {
elm[d:ndose] = 1
if (d == 1) {
earlystop = 1
break
}
}
if (extrasafe) {
if (d == 1 && n[1] >= 3) {
if (1 - pbeta(target, y[1] + 0.05, n[1] -
y[1] + 0.1) > cutoff.eli - offset) {
earlystop = 1
break
}
}
}
}
if (y[d] <= b.e[n[d]] && d != ndose) {
if (elm[d + 1] == 0)
d = d + 1
}
else if (y[d] >= b.d[n[d]] && d != 1) {
d = d - 1
}
else {
d = d
}
if (n[d] >= n.earlystop)
break
if (sum(n) >= (ncohort * cohortsize))
break
}
if (earlystop == 1) {
selectdose = 99
elm = rep(1, ndose)
}
else {
wsmtd = waterfall.subtrial.mtd(target, n, y, cutoff.eli,
extrasafe, offset, boundMTD=boundMTD, p.tox,temp)
selectdose = wsmtd$selectdose
is.escalation = wsmtd$is.escalation
}
npts[dosespace] = n
ntox[dosespace] = y
elimi[dosespace] = elm
list(ncohort = icohort, ntotal = icohort * cohortsize,
startdose = startdose, npts = npts, ntox = ntox,
totaltox = sum(ntox), totaln = sum(npts), pctearlystop = sum(selectdose ==
99) * 100, selectdose = selectdose, is.escalation = is.escalation,
elimi = elimi)
}
get.oc.comb.waterfall <- function(p.true, target, ncohort,
cohortsize, n.earlystop = 12, cutoff.eli = 0.95, p.saf = 0.6 *
target, p.tox = 1.4 * target, titration = FALSE,
extrasafe = FALSE, offset = 0.05, boundMTD=FALSE,ntrial = 1000) {
temp = get.boundary(target, ncohort = 150, cohortsize = 1,
cutoff.eli = cutoff.eli, extrasafe = extrasafe)$boundary_tab
JJ = nrow(p.true)
KK = ncol(p.true)
if (JJ > KK)
p.true = t(p.true)
true.mtd.position = cbind(1:JJ, apply(p.true, 1, function(x) {
flaga = rep(0, length(x))
tmp = which.min(abs(x - target))
flaga[tmp] = 1
flagb = (x <= target + 0.05)
ifelse(sum(flaga & flagb) > 0, which(flaga & flagb),
99)
}))
true.mtd.pos = apply(true.mtd.position, 1, function(x) (x[2] -
1) * JJ + x[1])
true.mtd.pos.new = true.mtd.pos[true.mtd.pos <= JJ *
KK]
nMTDs = paste(sort(true.mtd.pos.new), collapse = ",")
greater.than.contour = rep(1, length(JJ * KK))
less.than.contour = 1 - greater.than.contour
if (sum(true.mtd.position[, 2] <= KK) > 0) {
tmp = true.mtd.position[true.mtd.position[, 2] <=
JJ * KK, ]
if (length(tmp) == 2)
tmp = matrix(tmp, ncol = 2)
ONES = matrix(1, nrow = JJ, ncol = KK)
for (kkk in 1:nrow(tmp)) {
tmpa = tmp[kkk, 1]
tmpb = tmp[kkk, 2]
ONES[1:tmpa, 1:tmpb] = 0
}
ONES[true.mtd.pos.new] = 1
less.than.contour = which(ONES == 0)
ONES = matrix(1, nrow = JJ, ncol = KK)
for (kkk in 1:nrow(tmp)) {
tmpa = tmp[kkk, 1]
tmpb = tmp[kkk, 2]
ONES[tmpa:JJ, tmpb:KK] = 0
}
ONES[true.mtd.pos.new] = 1
greater.than.contour = which(ONES == 0)
}
at.contour = c(1:(JJ * KK))[-c(greater.than.contour,
less.than.contour)]
aa = function(x) as.numeric(as.character(x))
if (target < 0.05) {
stop(" the target is too low! ")
}
if (target > 0.6) {
stop("the target is too high!")
}
if ((target - p.saf) < (0.1 * target)) {
stop("the probability deemed safe cannot be higher than or too close to the target! ")
}
if ((p.tox - target) < (0.1 * target)) {
stop("the probability deemed toxic cannot be lower than or too close to the target!")
}
ndoses1 <- nrow(p.true)
ndoses2 <- ncol(p.true)
ntrial.phase1 = NULL
ntrial.mtd = NULL
ntrial.nt = NULL
ntrial.yt = NULL
ntrial.ne = NULL
ntrial.ye = NULL
if (cohortsize == 1)
titration = FALSE
titration.first.trial = FALSE
for (trial in 1:ntrial) {
trial.result = NULL
ntox = matrix(rep(0, (ndoses2) * ndoses1), ncol = ndoses2)
colnames(ntox) = paste("ntoxDoseB", 1:ndoses2, sep = "")
npts = matrix(rep(0, (ndoses2) * ndoses1), ncol = ndoses2)
colnames(npts) = paste("nptsDoseB", 1:ndoses2, sep = "")
elimi = matrix(0, nrow = ndoses1, ncol = ndoses2)
colnames(elimi) = paste("elimiDoseB", 1:ndoses2,
sep = "")
mtd = cbind(selectdoseA = 1:ndoses1, selectdoseB = rep(NA,
ndoses1))
trial.result = data.frame(cbind(trial = rep(trial,
ndoses1), mtd, npts, ntox, elimi))
totaln = 0
startdose = 1
dosespace = c(1:(ndoses1 - 1), (1:ndoses2) * ndoses1)
subtriali = 1
while (totaln < sum(ncohort) * cohortsize) {
if (titration & subtriali == 1) {
titration.first.trial = TRUE
}
else {
titration.first.trial = FALSE
}
subtrial = waterfall.subtrial(target, p.true = p.true,
dosespace = dosespace, npts = npts, ntox = ntox,
elimi = elimi, ncohort = ncohort[subtriali],
cohortsize = cohortsize, n.earlystop = n.earlystop,
startdose = startdose, p.saf = p.saf, p.tox = p.tox,
cutoff.eli = cutoff.eli, extrasafe = extrasafe,
offset = offset, totaln = totaln,
titration.first.trial = titration.first.trial,
temp = temp,boundMTD=boundMTD)
selectdose = ifelse(subtrial$selectdose == 99,
99, dosespace[subtrial$selectdose])
if (selectdose == 99)
break
dj = ifelse(selectdose%%ndoses1 == 0, selectdose%/%ndoses1,
selectdose%/%ndoses1 + 1)
di = selectdose - (dj - 1) * ndoses1
totaln = aa(subtrial$totaln)
npts = subtrial$npts
ntox = subtrial$ntox
elimi = subtrial$elimi
if ((subtriali == 1) & (selectdose < ndoses1)) {
for (a in (di + 1):ndoses1) for (b in 1:ndoses2) elimi[a,
b] = 1
if (subtrial$is.escalation == 1) {
startdose = 1
dosespace1 = c(di + ((2:ndoses2) - 1) * ndoses1)
subtriali = subtriali + 1
subtrial1 = waterfall.subtrial(target, p.true = p.true,
dosespace = dosespace1, npts = npts, ntox = ntox,
elimi = elimi, ncohort = ncohort[subtriali],
cohortsize = cohortsize, n.earlystop = n.earlystop,
startdose = startdose, p.saf = p.saf, p.tox = p.tox,
cutoff.eli, extrasafe, offset, totaln = totaln,
temp = temp,boundMTD=boundMTD)
selectdose1 = ifelse(subtrial1$selectdose ==
99, selectdose, dosespace1[subtrial1$selectdose])
if (selectdose1 == 99)
break
dj = ifelse(selectdose1%%ndoses1 == 0, selectdose1%/%ndoses1,
selectdose1%/%ndoses1 + 1)
di = selectdose1 - (dj - 1) * ndoses1
totaln = aa(subtrial1$totaln)
npts = subtrial1$npts
ntox = subtrial1$ntox
elimi = subtrial1$elimi
}
}
if (di - 1 == 0)
break
subtriali = subtriali + 1
if (dj < ndoses2)
elimi[di, (dj + 1):ndoses2] = 1
startdose = dj
dosespace = di - 1 + ((2:ndoses2) - 1) * ndoses1
if (dj == ndoses2)
startdose = dj - 1
}
npts = t(apply(npts, 1, aa))
ntox = t(apply(ntox, 1, aa))
elimi = t(apply(elimi, 1, aa))
phat = (ntox + 0.05)/(npts + 0.1)
phat = t(apply(phat, 1, aa))
colnames(phat) = paste("phat", 1:ndoses2, sep = "")
phat[elimi == 1] = 1.1
phat = Iso::biviso(phat, npts + 0.1, warn = TRUE)[,
]
phat = phat + (1e-05) * (matrix(rep(1:dim(npts)[1],
each = dim(npts)[2], len = length(npts)), dim(npts)[1],
byrow = T) + matrix(rep(1:dim(npts)[2], each = dim(npts)[1],
len = length(npts)), dim(npts)[1]))
colnames(phat) = paste("phat", 1:ndoses2, sep = "")
for (k in ndoses1:1) {
kn = npts[k, ]
ky = ntox[k, ]
kelimi = elimi[k, ]
kphat = phat[k, ]
if (kelimi[1] == 1 || sum(npts[kelimi == 0]) ==
0) {
kseldose = 99
}
else {
adm.set = (kn != 0) & (kelimi == 0)
adm.index = which(adm.set == T)
y.adm = ky[adm.set]
n.adm = kn[adm.set]
selectd = sort(abs(kphat[adm.set] - target),
index.return = T)$ix[1]
kseldose = adm.index[selectd]
}
mtd[k, 2] = kseldose
if (k < ndoses1)
if (mtd[k, 2] - mtd[k + 1, 2] <= 0 & mtd[k +
1, 2] != 99)
mtd[k, 2] = mtd[k + 1, 2]
}
trial.result[1:ndoses1, grep("nptsDoseB", colnames(trial.result))] = npts
trial.result[1:ndoses1, grep("ntoxDoseB", colnames(trial.result))] = ntox
trial.result[1:ndoses1, grep("selectdose", colnames(trial.result))] = mtd
trial.result[1:ndoses1, grep("elimiDoseB", colnames(trial.result))] = elimi
ntrial.mtd = rbind(ntrial.mtd, cbind(trial = rep(trial,
nrow(mtd)), mtd))
ntrial.nt = rbind(ntrial.nt, cbind(trial = rep(trial,
nrow(npts)), npts))
ntrial.yt = rbind(ntrial.yt, cbind(trial = rep(trial,
nrow(ntox)), ntox))
trial.result = cbind(trial.result, phat)
ntrial.phase1 = rbind(ntrial.phase1, trial.result)
}
ntrial.mtd = data.frame(ntrial.mtd)
colnames(ntrial.mtd) = c("trial", "doseA", "doseB")
alltoxpercent = round(sum(sapply(1:ntrial, function(x) ifelse(sum(ntrial.mtd$doseB[ntrial.mtd$trial ==
x] == 99) == ndoses1, 1, 0)), na.rm = T) * 100/1000,
3)
stoppercent = round(sum(sapply(1:ntrial, function(x) sum(ntrial.mtd$doseB[ntrial.mtd$trial ==
x] == 99))) * 100/1000, 1)
selpercent = matrix(0, nrow = ndoses1, ncol = ndoses2)
nselpercent = 0
selpercent1 = 0
selpercent2 = 0
mtdtable = NULL
mtdlist = list()
for (triali in 1:ntrial) {
mtddata = unique(as.matrix(ntrial.mtd[ntrial.mtd$trial ==
triali, 2:3]))
mtddata = mtddata[!is.na(mtddata[, 2]) & mtddata[,
2] <= KK, ]
if (length(mtddata) > 0) {
if (length(mtddata) == 2)
mtddata = matrix(mtddata, ncol = 2)
mtdlevel = aa(t(apply(mtddata, 1, function(x) (x[2] -
1) * ndoses1 + x[1])))
mtdlevel[mtdlevel > ndoses1 * ndoses2] = 99
if (sum(mtdlevel <= ndoses1 * ndoses2) > 0) {
selpercent[mtdlevel[mtdlevel <= ndoses1 * ndoses2]] = selpercent[mtdlevel[mtdlevel <=
ndoses1 * ndoses2]] + 1
mtdlist[[triali]] = mtdlevel[mtdlevel <= ndoses1 *
ndoses2]
mtdtable = c(mtdtable, paste(sort(aa(mtdlevel[mtdlevel <=
ndoses1 * ndoses2])), collapse = ","))
if (paste(sort(aa(mtdlevel[mtdlevel <= ndoses1 *
ndoses2])), collapse = ",") == nMTDs)
nselpercent = nselpercent + 1
if (sum(mtdlevel <= ndoses1 * ndoses2) == 1 &
sum(is.element(mtdlevel[mtdlevel <= ndoses1 *
ndoses2], which(p.true == target)) == F) ==
0)
selpercent1 = selpercent1 + 1
if (sum(mtdlevel <= ndoses1 * ndoses2) == 2 &
sum(is.element(mtdlevel[mtdlevel <= ndoses1 *
ndoses2], which(p.true == target)) == F) ==
0)
selpercent2 = selpercent2 + 1
}
else {
mtdlist[[triali]] = 99
mtdtable = c(mtdtable, 99)
}
}
}
selpercent = round(selpercent * 100/ntrial, 2)
if (JJ > KK) {
nptsdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.nt), by = ndoses1)) nptsdose = nptsdose +
ntrial.nt[i + 0:(ndoses1 - 1), -1]
ntoxdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.yt), by = ndoses1)) ntoxdose = ntoxdose +
ntrial.yt[i + 0:(ndoses1 - 1), -1]
colnames(p.true) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(p.true) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(nptsdose) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(nptsdose) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
colnames(ntoxdose) = paste("DoseB", 1:dim(t(p.true))[1],
sep = "")
rownames(ntoxdose) = paste("DoseA", 1:dim(t(p.true))[2],
sep = "")
out = list(p.true = round(t(p.true), 2), selpercent = round(t(selpercent),
2), npatients = round(t(nptsdose)/ntrial, 2),
ntox = round(t(ntoxdose)/ntrial, 2), totaltox = round(sum(ntoxdose/ntrial),
1), totaln = round(sum(nptsdose/ntrial), 1))
out2 = list()
if (length(at.contour) > 0) {
out2 = list(npercent.contour = paste(round(100 *
sum(nptsdose[at.contour])/sum(nptsdose), 1),
"%", sep = ""), npercent.above.contour = paste(round(100 *
sum(nptsdose[greater.than.contour])/sum(nptsdose),
1), "%", sep = ""), npercent.below.contour = paste(round(100 *
sum(nptsdose[less.than.contour])/sum(nptsdose),
1), "%", sep = ""), pcs.contour = paste(round(100 *
nselpercent/ntrial, 1), "%", sep = ""), flowchart = FALSE)
}
outnew = c(out, out2)
class(outnew)<-"boin"
return(outnew)
}
else {
nptsdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.nt), by = ndoses1)) nptsdose = nptsdose +
ntrial.nt[i + 0:(ndoses1 - 1), -1]
ntoxdose = matrix(0, nrow = ndoses1, ncol = ndoses2)
for (i in seq(1, nrow(ntrial.yt), by = ndoses1)) ntoxdose = ntoxdose +
ntrial.yt[i + 0:(ndoses1 - 1), -1]
out2 = list()
if (length(at.contour) > 0) {
out2 = list(npercent.contour = paste(round(100 *
sum(nptsdose[at.contour])/sum(nptsdose), 1),
"%", sep = ""), npercent.above.contour = paste(round(100 *
sum(nptsdose[greater.than.contour])/sum(nptsdose),
1), "%", sep = ""), npercent.below.contour = paste(round(100 *
sum(nptsdose[less.than.contour])/sum(nptsdose),
1), "%", sep = ""), pcs.contour = paste(round(100 *
nselpercent/ntrial, 1), "%", sep = ""))
}
rownames(p.true) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(p.true) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(selpercent) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(selpercent) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(nptsdose) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(nptsdose) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
rownames(ntoxdose) = paste("DoseA", 1:dim(p.true)[1],
sep = "")
colnames(ntoxdose) = paste("DoseB", 1:dim(p.true)[2],
sep = "")
out = list(p.true = apply(formatC(p.true, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
selpercent = apply(formatC(selpercent, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
npatients = apply(formatC(nptsdose/ntrial, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
ntox = apply(formatC(ntoxdose/ntrial, digits = 2,
format = "f", width = 5), c(1, 2), as.numeric),
totaltox = as.numeric(formatC(sum(ntoxdose/ntrial),
digits = 1, format = "f")), totaln = as.numeric(formatC(sum(nptsdose/ntrial),
digits = 1, format = "f")), flowchart = FALSE)
outnew = c(out, out2)
class(outnew)<-"boin"
return(outnew)
}
}
set.seed(seed)
JJ = nrow(p.true)
KK = ncol(p.true)
if (JJ > KK) {
stop("p.true should be arranged in a way (i.e., rotated) such that the number of rows is less than or equal to the number of columns.")
}
if (mtd.contour == FALSE) {
if (is.null(n.earlystop) == TRUE)
n.earlystop = 100
if (n.earlystop <= 6) {
warning("the value of n.earlystop is too low to ensure good operating characteristics. Recommend n.earlystop = 9 to 18 ")
}
if (length(ncohort) > 1) {
warning("ncohort is the total number of cohorts for the trial. Please enter a scalar.")
}
if (((JJ * KK) <= 4) & (sum(ncohort) <= 6)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
if (((JJ * KK) > 4) & (sum(ncohort) <= 8)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
out=get.oc.comb.boin(target = target, p.true = p.true,
ncohort = sum(ncohort), cohortsize = cohortsize,
n.earlystop = n.earlystop, startdose = startdose,
titration = titration,
p.saf = p.saf, p.tox = p.tox,
cutoff.eli = cutoff.eli, extrasafe = extrasafe,
offset = offset,boundMTD=boundMTD,
ntrial = ntrial)
class(out)<-"boin"
return(out)
}
if (mtd.contour == TRUE) {
if (missing(ncohort) == TRUE) {
constSeq = rep(round(4/cohortsize, 2), 20)
dosespaceSeq = c(JJ + KK - 1, rep(KK - 1, JJ - 1))
ncohort = ceiling(constSeq[1:JJ] * dosespaceSeq)
}
if (length(ncohort) != JJ) {
stop("The vector length of ncohort doesn't match the number of subtrials (the number of dose matrix rows).\n Please enter the number of ncohorts of each subtrial.")
}
if (((JJ * KK) <= 4) & (sum(ncohort) <= 6)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
if (((JJ * KK) > 4) & (sum(ncohort) <= 8)) {
warning("the sample size is too small, which may lead to poor operating characteristics. Suggest to increase the number of cohort.")
}
if (is.null(n.earlystop) == TRUE)
n.earlystop = 12
if (n.earlystop <= 6) {
warning("the value of n.earlystop is too low to ensure good operating characteristics. Recommend n.earlystop = 9 to 18 ")
}
out=get.oc.comb.waterfall(p.true = p.true, target = target,
ncohort = ncohort, cohortsize = cohortsize, n.earlystop = n.earlystop,
cutoff.eli = cutoff.eli, p.saf = p.saf, p.tox = p.tox,
titration = titration,
extrasafe = extrasafe, offset = offset,boundMTD=boundMTD,
ntrial = ntrial)
class(out)<-"boin"
return(out)
}
}
Try the BOIN package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.