bayer/R_DK/crmpack_R_programs_amended/ZZ_Rules-methods_DK.R
In 0liver0815/onc-crmpack-test: Object-Oriented Implementation of CRM Designs
#####################################################################################
## Author: Daniel Sabanes Bove [sabanesd *a*t* roche *.* com],
## Wai Yin Yeung [ w *.* yeung1 *a*t* lancaster *.* ac *.* uk]
## Project: Object-oriented implementation of CRM designs
##
## Time-stamp: <[Rules-methods.R] by DSB Die 09/06/2015 21:29>
##
## Description:
## Encapsulate the rule functions in formal methods.
##
## History:
## 07/02/2014 file creation
## 10/07/2015 Adding more classes
###################################################################################
##' @include Model-methods.R
##' @include Samples-class.R
##' @include Rules-class.R
{}
## ============================================================
## --------------------------------------------------
## Find out what is the next best dose
## --------------------------------------------------
##' Find the next best dose
##'
##' Compute the recommended next best dose.
##'
##' This function outputs the next best dose recommendation based on the
##' corresponding rule \code{nextBest}, the posterior \code{samples} from the
##' \code{model} and the underlying \code{data}.
##'
##' @param nextBest The rule, an object of class \code{\linkS4class{NextBest}}
##' @param doselimit The maximum allowed next dose. If this is an empty (length
##' 0) vector, then no dose limit will be applied in the course of dose
##' recommendation calculation, and a corresponding warning is given.
##' @param samples the \code{\linkS4class{Samples}} object
##' @param model The model input, an object of class \code{\linkS4class{Model}}
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots possible additional arguments without method dispatch
##' @return a list with the next best dose (element \code{value})
##' on the grid defined in \code{data}, and a plot depicting this recommendation
##' (element \code{plot}). In case of multiple plots also an element \code{singlePlots}
##' is included which returns the list of single plots, which allows for further
##' customization of these. Also additional list elements describing the outcome
##' of the rule can be contained.
##'
##' @export
##' @keywords methods
setGeneric("nextBest",
def=
function(nextBest, doselimit, samples, model, data, ...){
## there should be no default method,
## therefore just forward to next method!
standardGeneric("nextBest")
},
valueClass="list")
## --------------------------------------------------
## The MTD method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the MTD rule
##'
##' @example examples/Rules-method-NextBestMTD.R
##' @importFrom ggplot2 ggplot geom_density xlab ylab xlim aes geom_vline
##' geom_text
setMethod("nextBest",
signature=
signature(nextBest="NextBestMTD",
doselimit="numeric",
samples="Samples",
model="Model",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## First, generate the MTD samples.
mtdSamples <- dose(prob=nextBest@target,
model,
samples)
## then derive the next best dose
mtdEstimate <- nextBest@derive(mtdSamples=mtdSamples)
## be sure which doses are ok with respect to maximum
## possible dose - if one was specified
dosesOK <-
if(length(doselimit))
which(data@doseGrid <= doselimit)
else
seq_along(data@doseGrid)
## but now, round to the next possible grid point
index <- which.min(abs(data@doseGrid[dosesOK] - mtdEstimate))
ret <- data@doseGrid[dosesOK][index]
## produce plot
plot1 <- ggplot() +
geom_density(data=
data.frame(x=mtdSamples),
aes(x=x),
fill = "grey50", colour = "grey50") +
xlab("MTD") + ylab("Posterior density") +
xlim(range(data@doseGrid))
plot1 <- plot1 +
geom_vline(xintercept=mtdEstimate, colour="black", lwd=1.1) +
geom_text(data=
data.frame(x=mtdEstimate),
aes(x, 0,
label = "Est", hjust=+1, vjust = +1),
colour="black")
if(length(doselimit))
{
plot1 <- plot1 +
geom_vline(xintercept=doselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=doselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = +1),
colour="red")
}
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0, vjust = +1),
colour="blue")
## return next best dose and plot
return(list(value=ret,
plot=plot1))
})
## --------------------------------------------------
## The NCRM method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the NCRM method. The
##' additional list element \code{probs} contains the target and overdosing
##' probabilities (across all doses in the dose grid)
##' used in the derivation of the next best dose.
##'
##' @example examples/Rules-method-NextBestNCRM.R
##' @importFrom ggplot2 ggplot geom_bar xlab ylab ylim aes geom_vline
##' geom_hline geom_point
##' @importFrom gridExtra arrangeGrob
setMethod("nextBest",
signature=
signature(nextBest="NextBestNCRM",
doselimit="numeric",
samples="Samples",
model="Model",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## first we have to get samples from the dose-tox
## curve at the dose grid points.
probSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the probs, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
probSamples[, i] <- prob(dose=data@doseGrid[i],
model,
samples)
}
## Now compute probabilities to be in target and
## overdose tox interval
probTarget <-
colMeans((probSamples >= nextBest@target[1]) &
(probSamples <= nextBest@target[2]))
probOverdose <-
colMeans((probSamples > nextBest@overdose[1]) &
(probSamples <= nextBest@overdose[2]))
## which doses are eligible after accounting
## for maximum possible dose and discarding overdoses?
dosesBelowLimit <-
if(length(doselimit))
(data@doseGrid <= doselimit)
else
rep(TRUE, length(data@doseGrid))
dosesOK <- which(dosesBelowLimit &
(probOverdose < nextBest@maxOverdoseProb))
## check if there are doses that are OK
if(length(dosesOK))
{
## what is the recommended dose level?
## if maximum target probability is higher than some numerical
## threshold, then take that level, otherwise stick to the
## maximum level that is OK:
doseLevel <-
if(max(probTarget[dosesOK]) > 0.05)
{
which.max(probTarget[dosesOK])
} else {
which.max(data@doseGrid[dosesOK])
}
ret <- data@doseGrid[dosesOK][doseLevel]
} else {
## if none of the doses is OK:
doseLevel <- NA
ret <- NA
}
## produce plot
## first for the target probability
plot1 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probTarget * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="darkgreen",
fill="darkgreen") +
xlab("Dose") +
ylab(paste("Target probability [%]")) +
ylim(c(0, 100))
if(length(doselimit))
{
plot1 <- plot1 +
geom_vline(xintercept=doselimit,
lwd=1.1,
lty=2,
colour="black")
}
if(length(dosesOK))
{
plot1 <- plot1 +
geom_vline(xintercept=data@doseGrid[max(dosesOK)],
lwd=1.1,
lty=2,
colour="red")
plot1 <- plot1 +
geom_point(data=
data.frame(x=ret,
y=probTarget[dosesOK][doseLevel] *
100 + 0.03),
aes(x=x, y=y),
size=3,
pch=25,
col="red",
bg="red")
}
## second for the overdosing probability
plot2 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probOverdose * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="red",
fill="red") +
xlab("Dose") +
ylab("Overdose probability [%]") +
ylim(c(0, 100))
plot2 <- plot2 +
geom_hline(yintercept=nextBest@maxOverdoseProb * 100,
lwd=1.1,
lty=2,
colour="black")
## now plot them below each other
plotJoint <- gridExtra::arrangeGrob(plot1, plot2, nrow=2)
## return value and plot
return(list(value=ret,
plot=plotJoint,
singlePlots=list(plot1=plot1,
plot2=plot2),
probs=cbind(dose=data@doseGrid,
target=probTarget,
overdose=probOverdose)))
})
##' @describeIn nextBest Find the next best dose based on the NCRM method when
##' two parts trial is used - todo: need an example here for DataParts
##' @example examples/Rules-method-NextBestNCRM-DataParts.R
setMethod("nextBest",
signature=
signature(nextBest="NextBestNCRM",
doselimit="numeric",
samples="Samples",
model="Model",
data="DataParts"),
def=
function(nextBest, doselimit, samples, model, data, ...){
## exception when we are in part I or about to start part II!
if(all(data@part == 1L))
{
## here we will always propose the highest possible dose
## (assuming that the dose limit came from reasonable
## increments rule, i.e. inrementsRelativeParts)
if(identical(length(doselimit), 0L))
{
stop("doselimit needs to be given for Part I")
}
return(list(value=doselimit,
plot=NULL))
} else {
## otherwise we will just do the standard thing
callNextMethod(nextBest, doselimit, samples, model, data, ...)
}
})
## --------------------------------------------------
## The 3+3 method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the 3+3 method
##' @example examples/Rules-method-NextBestThreePlusThree.R
setMethod("nextBest",
signature=
signature(nextBest="NextBestThreePlusThree",
doselimit="missing",
samples="missing",
model="missing",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
## split the DLTs into the dose level groups
dltSplit <- split(data@y,
factor(data@x,
levels=data@doseGrid))
## number of patients and number of DLTs per dose level group
nPatients <- sapply(dltSplit, length)
nDLTs <- sapply(dltSplit, sum)
## what was the last dose level tested?
lastLevel <- tail(data@xLevel, 1)
## if there are less than 1/3 DLTs at that level
if(nDLTs[lastLevel]/nPatients[lastLevel] < 1/3)
{
## we could escalate, unless this is the highest
## level or the higher level was tried already
## (in which case it was not safe)
if((lastLevel == length(data@doseGrid)) ||
(nPatients[lastLevel+1] > 0))
{
nextLevel <- lastLevel
} else {
nextLevel <- lastLevel + 1
}
} else if(nDLTs[lastLevel]/nPatients[lastLevel] > 1/3) {
## rate here is too high, therefore deescalate
nextLevel <- lastLevel - 1
} else {
## otherwise: rate is 1/3 == 2/6,
## then it depends on the number of patients:
## if more than 3, then deescalate, otherwise stay.
nextLevel <-
if(nPatients[lastLevel] > 3)
lastLevel - 1
else
lastLevel
}
## do we stop here? only if we have no MTD
## or the next level has been tried enough (more than
## three patients already)
stopHere <-
if(nextLevel == 0)
{
TRUE
} else {
nPatients[nextLevel] > 3
}
## return value and plot
return(list(value=
if(nextLevel == 0) NA else data@doseGrid[nextLevel],
stopHere=stopHere))
})
## --------------------------------------------------
## The method for the dual endpoint model
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the dual endpoint
##' model. The additional list element \code{probs} contains the target and
##' overdosing probabilities (across all doses in the dose grid) used in the
##' derivation of the next best dose.
##' @example examples/Rules-method-NextBestDualEndpoint.R
##' @importFrom ggplot2 ggplot geom_bar xlab ylab ylim aes geom_vline
##' geom_hline geom_point
##' @importFrom gridExtra arrangeGrob
setMethod("nextBest",
signature=
signature(nextBest="NextBestDualEndpoint",
doselimit="numeric",
samples="Samples",
model="DualEndpoint",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## get the biomarker level samples
## at the dose grid points.
biomLevelSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the biomLevels, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
biomLevelSamples[, i] <- biomLevel(dose=data@doseGrid[i],
xLevel=i,
model,
samples)
}
## biomLevelSamples <- samples@data$betaW
## now get samples from the dose-tox
## curve at the dose grid points.
probSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the probs, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
probSamples[, i] <- prob(dose=data@doseGrid[i],
model,
samples)
}
## if target is relative to maximum
if(nextBest@scale == "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@data)){
## For each sample, look which dose is maximizing the
## simultaneous probability to be in the target biomarker
## range and below overdose toxicity
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1,ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= nextBest@target[1] * samples@data$Emax &
biomLevelSamples[, x] <= nextBest@target[2] * samples@data$Emax) /
nrow(biomLevelSamples)
})
}else{
## For each sample, look which was the minimum dose giving
## relative target level
targetIndex <- apply(biomLevelSamples, 1L,
function(x){
rnx <- range(x)
min(which((x >= nextBest@target[1] * diff(rnx) + rnx[1]) &
(x <= nextBest@target[2] * diff(rnx) + rnx[1] + 1e-10))
)
})
probTarget <- numeric(ncol(biomLevelSamples))
tab <- table(targetIndex)
probTarget[as.numeric(names(tab))] <- tab
probTarget <- probTarget / nrow(biomLevelSamples)
}
} else {
## otherwise 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
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1, ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= nextBest@target[1] &
biomLevelSamples[, x] <= nextBest@target[2]) /
nrow(biomLevelSamples)
})
}
## Now compute probabilities to be in
## overdose tox interval
probOverdose <-
colMeans((probSamples > nextBest@overdose[1]) &
(probSamples <= nextBest@overdose[2]))
## which doses are eligible after accounting
## for maximum possible dose and discarding overdoses?
dosesBelowLimit <-
if(length(doselimit))
(data@doseGrid <= doselimit)
else
rep(TRUE, length(data@doseGrid))
dosesOK <- which(dosesBelowLimit &
(probOverdose < nextBest@maxOverdoseProb))
## check if there are doses that are OK
if(length(dosesOK))
{
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
## what is the recommended dose level?
## if maximum target probability is higher than the numerical
## threshold, then take that level, otherwise stick to the
## maximum level that is OK:
doseLevel <-
if(max(probTarget[dosesOK]) > nextBest@targetThresh)
{
which.max(probTarget[dosesOK])
} else {
which.max(data@doseGrid[dosesOK])
}
ret <- data@doseGrid[dosesOK][doseLevel]
} else {
## if none of the doses is OK:
doseLevel <- NA
ret <- NA
}
## produce plot
## first for the target probability
plot1 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probTarget * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="darkgreen",
fill="darkgreen") +
xlab("Dose") +
ylab(paste("Target probability [%]")) +
ylim(c(0, 100))
if(length(doselimit))
{
plot1 <- plot1 +
geom_vline(xintercept=doselimit,
lwd=1.1,
lty=2,
colour="black")
}
if(length(dosesOK))
{
plot1 <- plot1 +
geom_vline(xintercept=data@doseGrid[max(dosesOK)],
lwd=1.1,
lty=2,
colour="red")
plot1 <- plot1 +
geom_point(data=
data.frame(x=ret,
y=probTarget[dosesOK][doseLevel] *
100 + 0.03),
aes(x=x, y=y),
size=3,
pch=25,
col="red",
bg="red")
}
## second for the overdosing probability
plot2 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probOverdose * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="red",
fill="red") +
xlab("Dose") +
ylab("Overdose probability [%]") +
ylim(c(0, 100))
plot2 <- plot2 +
geom_hline(yintercept=nextBest@maxOverdoseProb * 100,
lwd=1.1,
lty=2,
colour="black")
## now plot them below each other
plotJoint <- gridExtra::arrangeGrob(plot1, plot2, nrow=2)
## return value and plot
return(list(value=ret,
plot=plotJoint,
singlePlots=list(plot1=plot1,
plot2=plot2),
probs=cbind(dose=data@doseGrid,
target=probTarget,
overdose=probOverdose)))
})
## --------------------------------------------------
## The CRM MTD method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on CRM MTD allocation
##' probability
##'
##' @importFrom ggplot2 ggplot geom_density xlab ylab xlim aes geom_vline
##' geom_text
setMethod("nextBest",
signature=
signature(nextBest="NextBestMTDCRM",
doselimit="numeric",
samples="Samples",
model="Model",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## first we have to get samples from the dose-tox
## curve at the dose grid points.
probSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the probs, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
probSamples[, i] <- prob(dose=data@doseGrid[i],
model,
samples)
}
## count the total number of samples
nrow <- sampleSize(samples@options)
## Now identify dose rate probabilities that are less than or
## equal to target
probTarget <-(probSamples <= nextBest@target)
## sum the number of doses for each row, that meet the
## criterion (TRUE)
sumMTD <- rowSums(probTarget[]=='TRUE')
## count the number of different values of the sums
countsumMTD <- table(sumMTD)
## if the dimension of countsumMTD ias not the same as the nGrid+1
## we need to impute with zero's the missing dose levels
if (dim(countsumMTD)<(data@nGrid+1)){
frame<-as.data.frame(countsumMTD)
framefull<-data.frame(sumMTD=0:data@nGrid)
frame1 <- merge(x=framefull ,y= frame,by="sumMTD", all.x = TRUE)
frame1$Freq[is.na(frame1$Freq)] = 0
countsumMTD <- as.table(frame1$Freq)
names(countsumMTD)<-c(0:data@nGrid)
}
## percentage of different values of the sums
freqMTD <- (countsumMTD/nrow)*100
## when the data object definition does NOT contain a placebo dose
if(data@placebo == "FALSE")
{
## if raw frequencies corresponding to planned doses are used for the nexbest dose
if(nextBest@pbomethod == "none")
{
countsumMTDnone<-countsumMTD[2:dim(countsumMTD)]
rr <- (countsumMTDnone/sum(countsumMTDnone))*100
rrmaxvalue<-rr[rr[]==max(rr)]
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## if cumulative frequencies corresponding to planned doses are used for the nexbest dose
## and the frequency of zero dose is added to the frequency of the lowest planned dose
if(nextBest@pbomethod == "min")
{
freqMTDdmin <-freqMTD
freqMTDdmin[2]<-freqMTDdmin[1]+freqMTDdmin[2]
rr<-freqMTDdmin[2:dim(freqMTDdmin)]
rrmaxvalue<-rr[rr[]==max(rr)]
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## if cumulative frequencies corresponding to planned doses are used for the nexbest dose
## and the frequency of zero dose is added to the frequency of the highest planned dose
if(nextBest@pbomethod == "max")
{
freqMTDdmax <-freqMTD
freqMTDdmax[dim(freqMTDdmax)]<-freqMTDdmax[1]+freqMTDdmax[dim(freqMTDdmax)]
rr<-freqMTDdmax[2:dim(freqMTDdmax)]
rrmaxvalue<-rr[rr[]==max(rr)]
}
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## when the data object definition contains a placebo dose
## if TRUE the first dose level in the grid is considered as PLACEBO
if(data@placebo == "TRUE")
{
rr<-freqMTD
rrmaxvalue<-rr[rr[]==max(rr)]
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## be sure which doses are ok with respect to maximum
## possible dose - if one was specified
dosesOK <-
if(length(doselimit)){
which(data@doseGrid <= doselimit)
}else{
seq_along(data@doseGrid)
}
dosesOKmax<-max(data@doseGrid[dosesOK])
if (doseMTD >= dosesOKmax){
ret <- dosesOKmax
}
if (doseMTD < dosesOKmax) {
ret <- doseMTD
}
## produce plot
plot1 <- ggplot() +
geom_line(data=data.frame(x=data@doseGrid,
y=as.numeric(rr) ),
aes(x=x, y=y),
linetype = "dashed",
size=1,
colour="red"
) +
geom_point(data=data.frame(x=data@doseGrid,
y=as.numeric(rr)
),
aes(x=x, y=y),
size=2,
colour="red"
) +
geom_text(data=data.frame(x=data@doseGrid,
y=as.numeric(rr)
),
mapping=aes(x=data@doseGrid,
y=as.numeric(rr) ,
label=paste(round(as.numeric(rr), 1),"%")
),
size = 3,
hjust = 0.5,
vjust = -1
)+
geom_vline(xintercept = doseMTD,
linetype = "dashed",
size=1,
colour="blue"
) +
annotate(geom = "text",
x = (doseMTD),
y = 97.5,
label = "Max",
colour="blue"
) +
geom_vline(xintercept = ret,
linetype = "dashed",
size=1,
colour="green") +
annotate(geom = "text",
x = ret,
y = 92.5,
colour="green",
label = "Next"
) +
xlab("Dose") +
ylab(paste("Allocation criterion [%]")) +
scale_x_continuous(limits=range(data@doseGrid) , breaks=data@doseGrid) +
scale_y_continuous(limits=c(0, 100) , breaks=seq(0,100,10))
## if doselimit exists add the dose range on the plot
if(length(doselimit))
{
plot1 <- plot1+
geom_rect(aes(xmin = data@doseGrid[1],
ymin = 0,
xmax =doselimit,
ymax = 100
),
alpha = 0.5
) +
annotate(geom = "text",
x = doselimit,
y = 95,
label = "Doselimit"
)
}
## return next best dose
return(list(value=ret,
plot=plot1))
})
## ============================================================
## --------------------------------------------------
## Determine the maximum possible next dose
## --------------------------------------------------
##' Determine the maximum possible next dose
##'
##' Determine the upper limit of the next dose based on the increments rule.
##'
##' This function outputs the maximum possible next dose, based on the
##' corresponding rule \code{increments} and the \code{data}.
##'
##' @param increments The rule, an object of class
##' \code{\linkS4class{Increments}}
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots further arguments
##' @return the maximum possible next dose
##'
##' @export
##' @keywords methods
setGeneric("maxDose",
def=
function(increments, data, ...){
## there should be no default method,
## therefore just forward to next method!
standardGeneric("maxDose")
},
valueClass="numeric")
## --------------------------------------------------
## The maximum allowable relative increments in intervals method
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' relative increments
##'
##' @example examples/Rules-method-maxDose-IncrementsRelative.R
setMethod("maxDose",
signature=
signature(increments="IncrementsRelative",
data="Data"),
def=
function(increments, data, ...){
## determine what was the last dose
lastDose <- tail(data@x, 1)
## determine in which interval this dose was
lastInterval <-
findInterval(x=lastDose,
vec=increments@intervals)
## so the maximum next dose is
ret <-
(1 + increments@increments[lastInterval]) *
lastDose
return(ret)
})
## --------------------------------------------------
## The maximum allowable number of dose levels method
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' maximum dose levels to increment for the next dose
##'
##' @example examples/Rules-method-maxDose-IncrementsNumDoseLevels.R
setMethod("maxDose",
signature=
signature(increments="IncrementsNumDoseLevels",
data="Data"),
def=
function(increments, data, ...){
## determine what was the level of the last dose
lastDoseLevel <- tail(data@xLevel, 1)
## determine the maximum next dose level
maxNextDoseLevel <- min(length(data@doseGrid),
lastDoseLevel + increments@maxLevels)
## so the maximum next dose is
ret <- data@doseGrid[maxNextDoseLevel]
return(ret)
})
## --------------------------------------------------
## The maximum allowable relative increments, with special rules for
## part 1 and beginning of part 2, method method
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' relative increments and part 1 and 2
##' @example examples/Rules-method-maxDose-IncrementsRelativeParts.R
setMethod("maxDose",
signature=
signature(increments="IncrementsRelativeParts",
data="DataParts"),
def=
function(increments, data, ...){
## determine if there are already cohorts
## belonging to part 2:
alreadyInPart2 <- any(data@part == 2L)
## if so, we just call the next higher method
if(alreadyInPart2)
{
callNextMethod(increments, data, ...)
} else {
## otherwise we have special rules.
## what dose level (index) has the highest dose
## so far?
lastDoseLevel <- matchTolerance(max(data@x),
data@part1Ladder)
## determine the next maximum dose
ret <-
if(data@nextPart == 1L)
{
## here the next cohort will still be in part 1.
## Therefore we just make one step on the part 1 ladder:
data@part1Ladder[lastDoseLevel + 1L]
} else {
## the next cohort will start part 2.
## if there was a DLT so far:
if(any(data@y == 1L))
{
data@part1Ladder[lastDoseLevel + increments@dltStart]
} else {
## otherwise
if(increments@cleanStart > 0)
{
## if we want to start part 2 higher than
## the last part 1 dose, use usual increments
callNextMethod(increments, data, ...)
} else {
## otherwise
data@part1Ladder[lastDoseLevel + increments@cleanStart]
}
}
}
return(ret)
}
})
## --------------------------------------------------
## The maximum allowable relative increments in terms of DLTs
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' relative increments determined by DLTs so far
##'
##' @example examples/Rules-method-maxDose-IncrementsRelativeDLT.R
setMethod("maxDose",
signature=
signature(increments="IncrementsRelativeDLT",
data="Data"),
def=
function(increments, data, ...){
## determine what was the last dose
lastDose <- tail(data@x, 1)
## determine how many DLTs have occurred so far
dltHappened <- sum(data@y)
## determine in which interval this is
interval <-
findInterval(x=dltHappened,
vec=increments@DLTintervals)
## so the maximum next dose is
ret <-
(1 + increments@increments[interval]) *
lastDose
return(ret)
})
## --------------------------------------------------
## The maximum allowable relative increments in terms of DLTs
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' multiple increment rules (taking the minimum across individual increments).
##'
##' @example examples/Rules-method-maxDose-IncrementMin.R
setMethod("maxDose",
signature=
signature(increments="IncrementMin",
data="Data"),
def=
function(increments, data, ...){
## apply the multiple increment rules
individualResults <-
sapply(increments@IncrementsList,
maxDose,
data=data,
...)
## so the maximum increment is the minimum across the individual increments
ret <- min(individualResults)
return(ret)
})
## --------------------------------------------------
## The maximum allowable number of dose levels method based on Hard Safety Rule - Beta(a,b)
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' maximum eligible dose level to increment for the next dose
##' below the target toxicity, by a % probability based on the Beta(a,b) distribution.
##' this can be calculated for any number of patients in a cohort with a defined probability threshold of 90%
##' i.e. 1- pbeta(target, x+a, n-x+b) < 0.90
##' where x=DLTs, n=patients for a dose and use Beta(a=1,b=1) (or probably a=b=0.5)
##' if a dose satisfies the above criterio then eligible doses, are only the ones that that lower than that dose
##'
setMethod("maxDose",
signature=
signature(increments="IncrementsHSRBeta",
data="Data"),
def=
function(increments, data, ...){
#Load data dose vector and find the number of total subjects dosed at each dose
n<-table(data@x)
#Load data events vector at each dose and find the number of DLTs
yd<-table(data@y,data@x)
#check if both DLTs and no DLTs are present
if (nrow(yd)==2){
#Number of DLTs
y1<-yd[2,]
}
#check if only DLTs or only no DLTs are present
if (nrow(yd)==1){
#Number of DLTs in each scenario:
#if only DLTs
if (dimnames(yd)[1]==1){y1=n}
#if only no DLTs
if (dimnames(yd)[1]==0){y1=0}
}
#Calculate the probability of DLTs at each dose tested
# i.e. 1- pbeta(target, y+a, n-y+b) < 0.90
# where y=DLTs, n=patients for a dose and use Beta(a,b)
probTarget<-(1- pbeta(increments@target, y1+increments@a, n-y1+increments@b))
# find toxic doses
toxic<-probTarget>=increments@prob
# find safe doses
safe<-probTarget<increments@prob
# are there any safe doses
safegrid<-sum(safe[]=="TRUE")
# are there any toxic doses
toxicgrid<-sum(toxic[]=="TRUE")
# if only safe doses exist (toxic doses do not exist all) then all grid doses are eligible and max next dose is the maximum planned dose grid
if (safegrid[]>0 && toxicgrid[]==0){
maxdose<-max(data@doseGrid[])
}else if (safegrid[]>0 && toxicgrid[]>0){
# if safe and toxic doses both exist, then max next dose is the maximum safe dose
maxdose<-data@doseGrid[safegrid]
}else if (safegrid[]==0 && toxicgrid[]>0){
# if safe doses do not exist and all doses toxic
maxdose<-min(data@doseGrid)
# no dose is eligible so we assign the lowest (first) dose. This will force the next recommended dose to be the lowest dose and will trigger the stopping rule stoppingruleHSFBeta.
# therefore when using this increment rule it is necessary that you also use the relevant stoppingruleHSFBeta with the same setting. In this
# way the absence of safe doses including the lowest one will trigger trial to stop as long as these 2 are set the same.
}
## so the maximum next dose is
ret <- maxdose
return(ret)
})
## -----------------------------------------------------------------------------------------
## The maximum allowable number of dose levels method relative to the max dose tested so far
## -----------------------------------------------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' maximum dose levels to increment for the next dose
##'
setMethod("maxDose",
signature=
signature(increments="IncrementsNumDoseLevelsMaxTested",
data="Data"),
def=
function(increments, data, ...){
## determine what was the level of the max dose tested so far
maxDoseLevel <- max(data@xLevel)
## determine the maximum next dose level
maxNextDoseLevel <- min(length(data@doseGrid),
maxDoseLevel + increments@maxLevels)
## so the maximum next dose is
ret <- data@doseGrid[maxNextDoseLevel]
return(ret)
})
## ============================================================
## --------------------------------------------------
## "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@stopList <- c(e1@stopList,
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@stopList <- c(e1,
e2@stopList)
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@stopList <- c(e1@stopList,
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@stopList <- c(e1,
e2@stopList)
return(e2)
})
## --------------------------------------------------
## Stop the trial?
## --------------------------------------------------
##' Stop the trial?
##'
##' This function returns whether to stop the trial.
##'
##' @param stopping The rule, an object of class
##' \code{\linkS4class{Stopping}}
##' @param dose the recommended next best dose
##' @param samples the \code{\linkS4class{Samples}} object
##' @param model The model input, an object of class \code{\linkS4class{Model}}
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots additional arguments
##'
##' @return logical value: \code{TRUE} if the trial can be stopped, \code{FALSE}
##' otherwise. It should have an attribute \code{message} which gives the reason
##' for the decision.
##'
##' @export
##' @example examples/Rules-method-CombiningStoppingRulesAndOr.R
##' @keywords methods
setGeneric("stopTrial",
def=
function(stopping, dose, samples, model, data, ...){
## if the recommended next dose is NA,
## stop in any case.
if(is.na(dose))
{
return(structure(TRUE,
message="Recommended next best dose is NA"))
} else if(data@placebo && dose == min(data@doseGrid))
{
return(structure(TRUE,
message="Recommended next best dose is placebo dose"))
}
## there should be no default method,
## therefore just forward to next method!
standardGeneric("stopTrial")
},
valueClass="logical")
## --------------------------------------------------
## 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@stopList,
stopTrial,
dose=dose,
model=model,
data=data,
...)
} else {
lapply(stopping@stopList,
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))
})
## --------------------------------------------------
## 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@stopList,
stopTrial,
dose=dose,
model=model,
data=data,
...)
} else {
lapply(stopping@stopList,
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))
})
## --------------------------------------------------
## 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@stopList,
stopTrial,
dose=dose,
model=model,
data=data,
...)
} else {
lapply(stopping@stopList,
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))
})
## --------------------------------------------------
## 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))
})
## -------------------------------------------------------------
## 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 <- 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))
})
## --------------------------------------------------
## 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))
})
## --------------------------------------------------
## 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))
})
## --------------------------------------------------
## Stopping based on probability of target tox interval
## --------------------------------------------------
##' @describeIn stopTrial Stop based on probability of target tox interval
##'
##' @example examples/Rules-method-stopTrial-StoppingTargetProb.R
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetProb",
dose="numeric",
samples="Samples",
model="Model",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## first we have to get samples from the dose-tox
## curve at the dose.
probSamples <- prob(dose=dose,
model,
samples)
## Now compute probability to be in target interval
probTarget <-
mean((probSamples >= stopping@target[1]) &
(probSamples <= stopping@target[2]))
## so can we stop?
doStop <- probTarget >= stopping@prob
## generate message
text <-
paste("Probability for target toxicity is",
round(probTarget * 100),
"% for dose",
dose,
"and thus",
ifelse(doStop, "above", "below"),
"the required",
round(stopping@prob * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## --------------------------------------------------
## 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="Model",
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(prob=stopping@target,
model,
samples)
## what is the absolute threshold?
absThresh <- stopping@thresh * dose
## what is the probability to be above this dose?
prob <- 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, "above", "below"),
"the required",
round(stopping@prob * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## --------------------------------------------------
## Stopping based on probability of targeting biomarker
## --------------------------------------------------
##' @describeIn stopTrial Stop based on probability of targeting biomarker
##'
##' @example examples/Rules-method-stopTrial-StoppingTargetBiomarker.R
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetBiomarker",
dose="numeric",
samples="Samples",
model="DualEndpoint",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## compute the target biomarker prob at this dose
## get the biomarker level samples
## at the dose grid points.
biomLevelSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the biomLevels, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
biomLevelSamples[, i] <- biomLevel(dose=data@doseGrid[i],
xLevel=i,
model,
samples)
}
## if target is relative to maximum
if(stopping@scale == "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@data)){
## For each sample, look which dose is maximizing the
## simultaneous probability to be in the target biomarker
## range and below overdose toxicity
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1,ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= stopping@target[1]*samples@data$Emax &
biomLevelSamples[, x] <= stopping@target[2]*samples@data$Emax) / nrow(biomLevelSamples)
})
}else{
## For each sample, look which was the minimum dose giving
## relative target level
targetIndex <- apply(biomLevelSamples, 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))
)
})
probTarget <- numeric(ncol(biomLevelSamples))
tab <- table(targetIndex)
probTarget[as.numeric(names(tab))] <- tab
probTarget <- probTarget / nrow(biomLevelSamples)
}
} else {
## otherwise 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
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1, ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= stopping@target[1] &
biomLevelSamples[, x] <= stopping@target[2]) /
nrow(biomLevelSamples)
})
}
## so for this dose we have:
probTarget <- probTarget[which(data@doseGrid == dose)]
## so can we stop?
doStop <- probTarget >= stopping@prob
## generate message
text <-
paste("Probability for target biomarker is",
round(probTarget * 100),
"% for dose",
dose,
"and thus",
ifelse(doStop, "above", "below"),
"the required",
round(stopping@prob * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## --------------------------------------------------
## 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 <- (dose == data@doseGrid[data@nGrid])
return(structure(isHighestDose,
message=
paste("Next best dose is", dose, "and thus",
ifelse(isHighestDose, "the",
"not the"),
"highest dose")))
})
## -------------------------------------------------------------
## Stopping based on the lowest dose meeting the hard safety criteria using Beta based DLT probability
## -------------------------------------------------------------
##' @describeIn stopTrial Stop based on the lowest dose meeting the hard
##' safety criteria using Beta based DLT probability
setMethod("stopTrial",
signature=
signature(stopping="StoppingLowestDoseHSRBeta",
dose="numeric",
samples="ANY",
model="ANY",
data="Data"),
def=
function(stopping, dose, samples, model, data, ...){
#Load data dose vector and find the number of total subjects dosed at this the lowest planned dose
dose <- data@doseGrid[1]
if (dose %in% data@x){
xd<-table(data@x==dose)
if (dim(xd)>0){
n<-xd["TRUE"]
#Load data events vector at this dose and find the number of DLTs
yd<-table(data@y,data@x==dose)
#check if both DLTs and no DLTs are present
if (nrow(yd)==2){
#Number of DLTs
y1<-yd[2,"TRUE"]
}
#check if only DLTs or only no DLTs are present
if (nrow(yd)==1){
#Number of DLTs in each scenario:
#if only DLTs
if (dimnames(yd)[1]==1){y1=xd}
#if only no DLTs
if (dimnames(yd)[1]==0){y1=0}
}
#Calculate the probability of DLTs
# i.e. 1- pbeta(target, y+a, n-y+b) < 0.90
# where y=DLTs, n=patients for a dose and use Beta(a,b)
probTarget<-(1- pbeta(stopping@target, y1+stopping@a, n-y1+stopping@b))
## so can we stop?
doStop <- probTarget >= stopping@prob
## generate message
text <-
paste("Probability that the lowest dose",dose,"being toxic based on Beta(",
stopping@a,",",stopping@b,") is",probTarget,
"and thus",
ifelse(doStop, "above", "below"),
"the prespecified risk of",
stopping@prob)
## return both
return(structure(doStop,
message=text))
}
}else{
probTarget<-NA
## so can we stop?
doStop <- is.numeric(probTarget)
## generate message
text <-
paste("Lowest dose",dose,"never tested before and thus probability the lowest dose being toxic based on Beta(",
stopping@a,",",stopping@b,") is",probTarget
)
## return both
return(structure(doStop,
message=text))
}
})
## --------------------------------------------------
## Stopping based on precision of MTD calculated as CV(MTD)
## --------------------------------------------------
##' @describeIn stopTrial Stop based on precision of MTD calculated as CV(MTD)
setMethod("stopTrial",
signature=
signature(stopping="StoppingMTDCV",
dose="numeric",
samples="Samples",
model="Model",
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(prob=stopping@target,
model,
samples)
## CV of MTD derived based on MTD samples
MTDCV <- (mad(MTDSamples)/median(MTDSamples))
## so can we stop?
doStop <- ((MTDCV <= stopping@thresh) & (MTDCV >= 0))
## generate message
text <-
paste("CV of MTD is",
round(MTDCV * 100),"% and thus",
ifelse(doStop, "below", "above"),
"the required precision threshold of",
round(stopping@thresh * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## -----------------------------------------------------------------------------------------------
## Stopping based on probability of target tox interval and number of patients near to the highest dose
## -----------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on probability of target tox interval and
##' number of patients near to the highest dose
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetProbPatientsNearHighestDose",
dose="numeric",
samples="Samples",
model="Model",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## first we have to get samples from the dose-tox
## curve at the dose.
dose <- max(data@doseGrid)
probSamples <- prob(dose=dose,
model,
samples)
## Now compute probability to be in target interval
probTarget <-
mean((probSamples >= stopping@target[1]) &
(probSamples <= stopping@target[2]))
## so can we stop due to toxicity interval?
doStopprobTarget <- probTarget >= stopping@prob
lower <- (100 - stopping@percentage) / 100 * dose
upper <- (100 + stopping@percentage) / 100 * dose
## how many patients lie there?
nPatients <- sum((data@x >= lower) & (data@x <= upper))
## so can we stop due to sufficient number of patients?
doStopnPatients <- nPatients >= stopping@nPatients
## so can we stop overall?
doStop <- all(c(doStopprobTarget,doStopnPatients)=="TRUE")
## generate message
text <-
paste("Probability for target toxicity is ",
round(probTarget * 100),
" % for highest dose ",
dose,
" and thus ",
ifelse(doStopprobTarget, "above", "below"),
" the required ",
round(stopping@prob * 100),
"%"," and ",
nPatients,
" patients lie within ",
stopping@percentage,
"% of the highest dose ",
dose,
". This ",
ifelse(doStopnPatients, "reached", "is below"),
" the required ",
stopping@nPatients,
" patients",
sep="")
## return both
return(structure(doStop,
message=text))
})
## -----------------------------------------------------------------------------------------------
## Stopping based on probability of target tox interval and number of patients near to the lowest dose
## -----------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on probability of target tox interval and
##' number of patients near to the lowest dose
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetProbPatientsNearLowestDose",
dose="numeric",
samples="Samples",
model="Model",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## first we have to get samples from the dose-tox
## curve at the dose.
dose <- data@doseGrid[1]
probSamples <- prob(dose=dose,
model,
samples)
## Now compute probability to be in target interval
probTarget <-
mean((probSamples >= stopping@target[1]) &
(probSamples <= stopping@target[2]))
## so can we stop due to toxicity interval?
doStopprobTarget <- probTarget >= stopping@prob
lower <- (100 - stopping@percentage) / 100 * dose
upper <- (100 + stopping@percentage) / 100 * dose
## how many patients lie there?
nPatients <- sum((data@x >= lower) & (data@x <= upper))
## so can we stop due to sufficient number of patients?
doStopnPatients <- nPatients >= stopping@nPatients
## so can we stop overall?
doStop <- all(c(doStopprobTarget,doStopnPatients)=="TRUE")
## generate message
text <-
paste("Probability for target toxicity is",
round(probTarget * 100),
" % for lowest dose ",
dose,
" and thus ",
ifelse(doStopprobTarget, "above", "below"),
" the required ",
round(stopping@prob * 100),
"%"," and ",
nPatients,
" patients lie within ",
stopping@percentage,
" % of the lowest dose ",
dose,
". This ",
ifelse(doStopnPatients, "reached", "is below"),
" the required ",
stopping@nPatients,
" patients",
sep="")
## return both
return(structure(doStop,
message=text))
})
## ============================================================
## --------------------------------------------------
## "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(...))
})
## --------------------------------------------------
## Determine the size of the next cohort
## --------------------------------------------------
##' Determine the size of the next cohort
##'
##' This function determines the size of the next cohort.
##'
##' @param cohortSize The rule, an object of class
##' \code{\linkS4class{CohortSize}}
##' @param dose the next dose
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots additional arguments
##'
##' @return the size as integer value
##'
##' @export
##' @keywords methods
setGeneric("size",
def=
function(cohortSize, dose, data, ...){
## if the recommended next dose is NA,
## don't check and return 0
if(is.na(dose))
{
return(0L)
}
## there should be no default method,
## therefore just forward to next method!
standardGeneric("size")
},
valueClass="integer")
## --------------------------------------------------
## The dose range method
## --------------------------------------------------
##' @describeIn size Determine the cohort size based on the range into which the
##' next dose falls into
##'
##' @example examples/Rules-method-size-CohortSizeRange.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeRange",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## determine in which interval the next dose is
interval <-
findInterval(x=dose,
vec=cohortSize@intervals)
## so the cohort size is
ret <- cohortSize@cohortSize[interval]
return(ret)
})
## --------------------------------------------------
## The DLT range method
## --------------------------------------------------
##' @describeIn size Determine the cohort size based on the number of DLTs so
##' far
##'
##' @example examples/Rules-method-size-CohortSizeDLT.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeDLT",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## determine how many DLTs have occurred so far
dltHappened <- sum(data@y)
## determine in which interval this is
interval <-
findInterval(x=dltHappened,
vec=cohortSize@DLTintervals)
## so the cohort size is
ret <- cohortSize@cohortSize[interval]
return(ret)
})
## --------------------------------------------------
## Size based on maximum of multiple cohort size rules
## --------------------------------------------------
##' @describeIn size Size based on maximum of multiple cohort size rules
##' @example examples/Rules-method-size-CohortSizeMax.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeMax",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## evaluate the individual cohort size rules
## in the list
individualResults <-
sapply(cohortSize@cohortSizeList,
size,
dose=dose,
data=data,
...)
## summarize to obtain overall result
overallResult <- max(individualResults)
return(overallResult)
})
## --------------------------------------------------
## Size based on minimum of multiple cohort size rules
## --------------------------------------------------
##' @describeIn size Size based on minimum of multiple cohort size rules
##' @example examples/Rules-method-size-CohortSizeMin.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeMin",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## evaluate the individual cohort size rules
## in the list
individualResults <-
sapply(cohortSize@cohortSizeList,
size,
dose=dose,
data=data,
...)
## summarize to obtain overall result
overallResult <- min(individualResults)
return(overallResult)
})
## --------------------------------------------------
## Constant cohort size
## --------------------------------------------------
##' @describeIn size Constant cohort size
##' @example examples/Rules-method-size-CohortSizeConst.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeConst",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
return(cohortSize@size)
})
## --------------------------------------------------
## Cohort size based on the parts
## --------------------------------------------------
##' @describeIn size Cohort size based on the parts
##' @example examples/Rules-method-size-CohortSizeParts.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeParts",
dose="ANY",
data="DataParts"),
def=
function(cohortSize, dose, data, ...){
return(cohortSize@sizes[data@nextPart])
})
## ================================================================
## The nextBest method based only on DLE responses with samples
## ================================================================
##' @describeIn nextBest Find the next best dose based on the 'NextBestTDsamples'
##' class rule. This a method based only on the DLE responses and for
##' \code{\linkS4class{LogisticIndepBeta}} model class object involving DLE samples
##'
##' @importFrom ggplot2 ggplot geom_density xlab ylab xlim aes geom_vline
##' geom_text
##'
##' @example examples/Rules-method-nextbest_TDsamples.R
##'
##' @export
##' @keywords methods
setMethod("nextBest",
signature=
signature(nextBest="NextBestTDsamples",
doselimit="numeric",
samples="Samples",
model="LogisticIndepBeta",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
## First, generate the TDtargetDuringTrial (TDtarget During a Trial) and
## TDtargetEndOfTrial (TDtarget at the EndOfTrial) samples.
TDtargetDuringTrialSamples <- dose(prob=nextBest@targetDuringTrial,
model,
samples)
TDtargetEndOfTrialSamples <- dose(prob=nextBest@targetEndOfTrial,
model,
samples)
## then derive the prior/posterior mean of the above two samples
TDtargetDuringTrialEstimate <- as.numeric(nextBest@derive(TDsamples=TDtargetDuringTrialSamples))
TDtargetEndOfTrialEstimate <- as.numeric(nextBest@derive(TDsamples=TDtargetEndOfTrialSamples))
## be sure which doses are ok with respect to maximum
## possible dose
dosesOK <- which(data@doseGrid <= doselimit)
##Find the index of next dose in the doseGrid
##next dose is the dose level closest below the TDtargetDuringTrialEstimate
index <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level (in doseGrid) closest below the TDtargetEndOfTrialEstimate
index1 <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret1 <- data@doseGrid[dosesOK][index1]
CITDEOT <- as.numeric(quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975)))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
## produce plot
plot1 <- ggplot() +
geom_density(data=
data.frame(x=TDtargetDuringTrialSamples),
aes(x=x),
fill = "grey50", colour = "grey50")
plot1 <- plot1 +
geom_density(data=
data.frame(x=TDtargetEndOfTrialSamples),
aes(x=x),
fill = "grey50", colour = "violet") +
xlab("TD") + ylab("Posterior density") +
xlim(range(data@doseGrid))
plot1 <- plot1+
geom_vline(xintercept=TDtargetDuringTrialEstimate, colour="orange", lwd=1.1) +
annotate("text",label=paste(paste("TD",nextBest@targetDuringTrial*100),"Estimate"),
x=TDtargetDuringTrialEstimate,y=0,hjust=-0.1, vjust = -20,size=5,colour="orange")
plot1 <- plot1+
geom_vline(xintercept=TDtargetEndOfTrialEstimate, colour="violet", lwd=1.1) +
annotate("text",label=paste(paste("TD",nextBest@targetEndOfTrial*100),"Estimate"),
x=TDtargetEndOfTrialEstimate,y=0,hjust=-0.1, vjust = -25,size=5,colour="violet")
if (doselimit > max(data@doseGrid)){maxdoselimit<-max(data@doseGrid)} else {maxdoselimit <-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -35),
colour="red")
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0.1, vjust = -30),
colour="blue")
## return next best dose and plot
return(list(nextdose=ret,
targetDuringTrial=nextBest@targetDuringTrial,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
targetEndOfTrial=nextBest@targetEndOfTrial,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=ret1,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
plot=plot1))
})
## -------------------------------------------------------------------------------
## The nextBest method based only on DLE responses without samples
## -----------------------------------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the 'NextBestTD'
##' class rule. This a method based only on the DLE responses and for
##' \code{\linkS4class{LogisticIndepBeta}} model class object without DLE samples
##'
##' @param SIM internal command to notify if this method is used within simulations. Default as FALSE
##' @importFrom ggplot2 ggplot xlab ylab xlim aes geom_vline
##' geom_text
##'
##' @example examples/Rules-method-nextbest_TD.R
##'
##' @export
##' @keywords methods
setMethod("nextBest",
signature=
signature(nextBest="NextBestTD",
doselimit="numeric",
samples="missing",
model="LogisticIndepBeta",
data="Data"),
def=
function(nextBest, doselimit, model, data, SIM=FALSE,...){
##Find the target prob During Trial
targetDuringTrial<- nextBest@targetDuringTrial
##Find the target prob End of Trial
targetEndOfTrial<- nextBest@targetEndOfTrial
mylabel<-targetDuringTrial*100
label2 <-targetEndOfTrial*100
## Find the TD30 Estimate and TD(target) Estimate
TDtargetEndOfTrialEstimate <- dose(prob=targetEndOfTrial,
model)
TDEfourdg<-signif(TDtargetEndOfTrialEstimate,digits=4)
TDtargetDuringTrialEstimate<-dose(prob=targetDuringTrial,
model)
TDDfourdg<-signif(TDtargetDuringTrialEstimate,digits=4)
probDLE=prob(dose=data@doseGrid,
model=model)
## be sure which doses are ok with respect to maximum
## possible dose
dosesOK <- which(data@doseGrid <= doselimit)
##Find the index of next dose in the doseGrid
##next dose is the dose level closest below the TDtargetEstimate
index <- suppressWarnings(max(which((TDDfourdg- data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level (in doseGrid) closest below the TD30Estimate
index <- suppressWarnings(max(which((TDEfourdg - data@doseGrid[dosesOK]) >= 0)))
retTDE <- data@doseGrid[dosesOK][index]
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(targetEndOfTrial/(1-targetEndOfTrial))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval
CITDEOT <- c()
CITDEOT[2] <- exp(log(TDtargetEndOfTrialEstimate) + 1.96* sqrt(varEta))
CITDEOT[1] <- exp(log(TDtargetEndOfTrialEstimate) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData <- data.frame(dose=data@doseGrid,
probDLE=prob(dose=data@doseGrid,
model=model))
##make the plot
gdata <- with(plotData,
data.frame(x=dose,
y=probDLE,
group=rep("Estimated DLE",each=nrow(plotData)),
Type=factor(rep("Estimated DLE",nrow(plotData)),levels="Estimated DLE")))
plot1 <- ggplot(data=gdata, aes(x=x,y=y), group=group) +
xlab("Dose Levels")+
ylab(paste("Probability of DLE")) + ylim(c(0,1)) + xlim(c(0,max(data@doseGrid))) +
geom_line(colour=I("red"), size=1.5)
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
if ((TDDfourdg < LowestDose)|(TDDfourdg > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("TD",targetDuringTrial*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plot1<- plot1}
} else {plot1 <- plot1+
geom_point(data=data.frame(x=TDtargetDuringTrialEstimate,y=targetDuringTrial),aes(x=x,y=y),colour="orange", shape=15, size=8) +
annotate("text",label=paste(paste("TD",mylabel),"Estimate"),x=TDtargetDuringTrialEstimate+1,y=targetDuringTrial-0.2,size=5,colour="orange")}
if ((TDEfourdg < LowestDose)|(TDEfourdg > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("TD",targetEndOfTrial*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {plot1 <- plot1+
geom_point(data=data.frame(x=TDtargetEndOfTrialEstimate,y=0.3),aes(x=x,y=y),colour="violet", shape=16, size=8) +
annotate("text",label=paste(paste("TD",label2),"Estimate"),x=TDtargetEndOfTrialEstimate+1,y=targetEndOfTrial-0.1,size=5,colour="violet")}
if (doselimit > max(data@doseGrid)) {maxdoselimit <- max(data@doseGrid)} else {maxdoselimit<-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="brown", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -30),
colour="brown")
plot1 <-plot1 +
geom_vline(xintercept=ret, colour="purple", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0, vjust = -30),
colour="purple")
## return next best dose and plot
return(list(nextdose=ret,
targetDuringTrial=targetDuringTrial,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
targetEndOfTrial=targetEndOfTrial,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialatdoseGrid=retTDE,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
plot=plot1))
})
## ------------------------------------------------------------------------------------
## the nextBest method based on DLE and efficacy responses without DLE and efficacy samples
## -------------------------------------------------------------------------- ----------
##' @describeIn nextBest for slots \code{nextBest},\code{doselimit}, \code{data} and \code{SIM}. This is
##' a function to find the next best dose based on the 'NextBestMaxGain'
##' class rule. This a method based on the DLE responses and efficacy responses without DLE and
##' efficacy samples.
##'
##' @param Effmodel the efficacy model of \code{\linkS4class{ModelEff}} class object
##'
##' @importFrom ggplot2 ggplot xlab ylab xlim aes geom_vline
##' geom_text
##'
##' @example examples/Rules-method-nextbest_MaxGain.R
##'
##' @importFrom ggplot2 scale_colour_manual
##' @export
##' @keywords methods
setMethod("nextBest",
signature=
signature(nextBest="NextBestMaxGain",
doselimit="numeric",
samples="missing",
model="ModelTox",
data="DataDual"),
def=
function(nextBest,doselimit,model,data,Effmodel,SIM=FALSE,...){
stopifnot(is(Effmodel, "ModelEff"))
DuringTrialtargetprob <- nextBest@DLEDuringTrialtarget
EndOfTrialtargetprob <- nextBest@DLEEndOfTrialtarget
## Find the TDtarget Estimate for During Trial and End of trial
TDtargetEndOfTrialEstimate <- dose(prob=EndOfTrialtargetprob,model=model)
TDtargetDuringTrialEstimate<-dose(prob=DuringTrialtargetprob,model=model)
##Get all prob of DLE at all dose levels
probDLE=prob(dose=data@doseGrid,
model=model)
##Define gain function
Gainfun<-function(DOSE){
-gain(DOSE,DLEmodel=model,Effmodel=Effmodel)
}
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
##Find the dose which gives the maximum gain
Gstar<-(optim(LowestDose,Gainfun, method = "L-BFGS-B", lower=LowestDose,upper=max(data@doseGrid))$par)
##Find the maximum gain value
MaxGain<--(optim(LowestDose,Gainfun,method = "L-BFGS-B", lower=LowestDose,upper=max(data@doseGrid))$value)
## be sure which doses are ok with respect to maximum
## possible dose
dosesOK <- which(data@doseGrid <= doselimit)
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
##FIND the next dose which is the minimum between TDtargetDuringTrial and Gstar
nextdose<-min(TDtargetDuringTrialEstimate,Gstar)
##Find the dose level in doseGrid closest below nextdose
index <- suppressWarnings(max(which((signif(nextdose,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level in doseGrid closest below TDtargetEndOfTrial
indexE <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retE <- data@doseGrid[indexE]
##Find the dose level in doseGrid closest below TDtargetDuringTrial
indexD <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retD <- data@doseGrid[indexD]
##Find the dose level in doseGrid closest below Gstar
Gstarindex <- suppressWarnings(max(which((signif(Gstar,digits=4) - data@doseGrid[dosesOK]) >= 0)))
Gstarret <- data@doseGrid[Gstarindex]
if (data@placebo){
logGstar <- log(Gstar+Effmodel@c)} 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
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 <- t(deltaG)%*%covBETA%*%deltaG
##Find the upper and lower limit of the 95% credibility interval of Gstar
CIGstar <-c()
CIGstar[2] <- exp(logGstar + 1.96* sqrt(varlogGstar))
CIGstar[1] <- exp(logGstar - 1.96* sqrt(varlogGstar))
##The ratio of the upper to the lower 95% credibility interval
ratioGstar <- as.numeric(CIGstar[2]/CIGstar[1])
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(EndOfTrialtargetprob/(1-EndOfTrialtargetprob))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval of
##TDtargetEndOfTrial
CITDEOT <- c()
CITDEOT[2] <- exp(log(TDtargetEndOfTrialEstimate) + 1.96* sqrt(varEta))
CITDEOT[1] <- exp(log(TDtargetEndOfTrialEstimate) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData<-data.frame(dose=rep(data@doseGrid,3),
values=c(prob(dose=data@doseGrid,
model=model),
ExpEff(dose=data@doseGrid,
model=Effmodel),
gain(dose=data@doseGrid,
DLEmodel=model,
Effmodel=Effmodel)))
gdata<-with(plotData,
data.frame(x=dose,
y=values,
group=c(rep("p(DLE)",length(data@doseGrid)),
rep("Expected Efficacy",length(data@doseGrid)),
rep("Gain",length(data@doseGrid))),
Type=factor("Estimate",levels="Estimate")
))
plot1 <- ggplot(data=gdata, aes(x=x,y=y))+geom_line(aes(group=group,color=group),size=1.5)+
ggplot2::scale_colour_manual(name="curves",values=c("blue","green3","red"))+
xlab("Dose Level")+ xlim(c(0,max(data@doseGrid)))+
ylab(paste("Values")) + ylim(c(min(gdata$y),max(gdata$y)))
if ((signif(TDtargetEndOfTrialEstimate,4) < LowestDose) |(signif(TDtargetEndOfTrialEstimate,4) > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",EndOfTrialtargetprob*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <-plot1 + geom_point(data=data.frame(x=TDtargetEndOfTrialEstimate,y=EndOfTrialtargetprob),aes(x=x,y=y),colour="violet", shape=16, size=8) +
annotate("text",label=paste(paste("TD",EndOfTrialtargetprob*100),"Estimate"),x=TDtargetEndOfTrialEstimate-3,y=0.2,size=5,colour="violet")}
if ((signif(Gstar,4) < LowestDose) |(signif(Gstar,4) > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste("Estimated Gstar=",paste(Gstar," not within dose Grid")))} else {plot1 <- plot1}
} else {plot1 <- plot1 +
geom_point(data=data.frame(x=Gstar,y=MaxGain),aes(x=x,y=y),colour="green3", shape=17, size=8) +
annotate("text",label="Max Gain Estimate",x=Gstar,y=MaxGain-0.1,size=5,colour="green3")}
mylabel=format(DuringTrialtargetprob,digits=2)
if ((signif(TDtargetDuringTrialEstimate,4) < LowestDose) |(signif(TDtargetDuringTrialEstimate,4) > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",DuringTrialtargetprob*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_point(data=data.frame(x=signif(TDtargetDuringTrialEstimate,4),y=DuringTrialtargetprob),aes(x=x,y=y),colour="orange", shape=15, size=8) +
annotate("text",label=paste(paste("TD",DuringTrialtargetprob*100),"Estimate"),x=TDtargetDuringTrialEstimate+25,
y=DuringTrialtargetprob+0.01,size=5,colour="orange")
}
if (doselimit > max(data@doseGrid)) {maxdoselimit <- max(data@doseGrid)} else {maxdoselimit<-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="brown", lwd=1.1) +
annotate("text",label="Max",x=maxdoselimit-2,y=max(gdata$y),size=5,colour="brown")
plot1 <-plot1 +
geom_vline(xintercept=ret, colour="purple", lwd=1.1) +
annotate("text",label="Next", x=ret+1, y=max(gdata$y)-0.05,size=5,color="purple")
## return next best dose and plot
return(list(nextdose=ret,
DLEDuringTrialtarget=DuringTrialtargetprob,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
TDtargetDuringTrialAtDoseGrid=retD,
DLEEndOfTrialtarget=EndOfTrialtargetprob,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=retE,
GstarEstimate=Gstar,
GstarAtDoseGrid=Gstarret,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
CIGstar=CIGstar,
ratioGstar=ratioGstar,
plot=plot1))
})
## =====================================================================================
## the nextBest method based on DLE and efficacy responses with DLE and efficacy samples
## -------------------------------------------------------------------------- ----------
##' @describeIn nextBest for slots \code{nextBest},\code{doselimit}, \code{data} and \code{SIM}. This is
##' a function to find the next best dose based on the 'NextBestMaxGainSamples'
##' class rule. This a method based on the DLE responses and efficacy responses with DLE and
##' efficacy samples. Effmodel must be of class \code{\linkS4class{Effloglog}} or
##' \code{\linkS4class{EffFlexi}}.
##' @param Effsamples the efficacy samples of \code{\linkS4class{Samples}} class object
##'
##' @importFrom ggplot2 ggplot geom_histogram xlab ylab xlim aes geom_vline
##' geom_text
##' @example examples/Rules-method-nextbest_MaxGainSamples.R
setMethod("nextBest",
signature=
signature(nextBest="NextBestMaxGainSamples",
doselimit="numeric",
samples="Samples",
model="ModelTox",
data="DataDual"),
def=
function(nextBest,doselimit,samples,model,data,Effmodel,Effsamples,SIM=FALSE,...){
if(is(Effmodel, "Effloglog"))
{
##first get the probDLE samples
points <- data@doseGrid
probDLESamples <- matrix(nrow=sampleSize(samples@options),
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:
probDLESamples[, i] <- prob(dose=points[i],
model=model,
samples=samples)
}
probDLE <- apply(probDLESamples,2,FUN=nextBest@TDderive)
DuringTrialtargetprob <- nextBest@DLEDuringTrialtarget
EndOfTrialtargetprob <- nextBest@DLEEndOfTrialtarget
## Find the TDtarget samples for During Trial and End of trial
TDtargetEndOfTrialSamples <- dose(prob=EndOfTrialtargetprob,
model=model,
samples=samples)
TDtargetDuringTrialSamples<-dose(prob=DuringTrialtargetprob,
model=model,
samples=samples)
## Find the TDtarget Estimate for During Trial and End of trial
TDtargetEndOfTrialEstimate <- as.numeric(nextBest@TDderive(TDtargetEndOfTrialSamples))
## Ensure the estimate is within dose range
#TDtargetEndOfTrialEstimate <- min(TDtargetEndOfTrialEstimate,max(data@doseGrid))
TDtargetDuringTrialEstimate<-as.numeric(nextBest@TDderive(TDtargetDuringTrialSamples))
## Ensure the estimate is within dose range
#TDtargetDuringTrialEstimate <- min(TDtargetDuringTrialEstimate,max(data@doseGrid))
##we have to get samples from the gain values at all dose levels
ExpEffSamples <- matrix(nrow=sampleSize(Effsamples@options),
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:
ExpEffSamples[, i] <- ExpEff(dose=points[i],
model=Effmodel,
samples=Effsamples)
}
ExpEff <- apply(ExpEffSamples,2,FUN=nextBest@Gstarderive)
GainSamples <- matrix(nrow=sampleSize(samples@options),
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],
DLEmodel=model,
DLEsamples=samples,
Effmodel=Effmodel,
Effsamples=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]
##Obtain the Gstar estimate which is the nth percentile of the Gstar samples
Gstar <- as.numeric(nextBest@Gstarderive(GstarSamples))
##Ensure the estimate is within dose range
#Gstar <- min(Gstar,max(data@doseGrid))
gainvalues <- apply(GainSamples,2,FUN=nextBest@Gstarderive)
dosesOK <- which(data@doseGrid <= doselimit)
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
##FIND the next dose which is the minimum between TDtargetDuringTrial and Gstar
nextdose<-min(TDtargetDuringTrialEstimate,Gstar)
##Find the dose level in doseGrid closest below nextdose
index <- suppressWarnings(max(which((signif(nextdose,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level in doseGrid closest below TDtargetEndOfTrial
indexE <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retE <- data@doseGrid[indexE]
##Find the dose level in doseGrid closest below TDtargetDuringTrial
indexD <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retD <- data@doseGrid[indexD]
##Find the dose level in doseGrid closest below Gstar
Gstarindex <- suppressWarnings(max(which((signif(Gstar,digits=4) - data@doseGrid[dosesOK]) >= 0)))
Gstarret <- data@doseGrid[Gstarindex]
##Find the 95% credibility interval of Gstar and its ratio of the upper to the lower limit
CIGstar <- as.numeric(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 <- as.numeric(quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975)))
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData<-data.frame(dose=rep(data@doseGrid,3),
values=c(probDLE,
ExpEff,
gainvalues))
gdata<-with(plotData,
data.frame(x=dose,
y=values,
group=c(rep("p(DLE)",length(data@doseGrid)),
rep("Expected Efficacy",length(data@doseGrid)),
rep("Gain",length(data@doseGrid))),
Type=factor("Estimate",levels="Estimate")
))
## produce plot
plot1 <- ggplot() +
geom_histogram(data=
data.frame(x=GstarSamples),
aes(x=x),
fill = "darkgreen", colour = "green3", binwidth = 25) +
xlab("Gstar")+ xlim(c(0,max(data@doseGrid)))+
ylab("Posterior density")
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
if (signif(TDtargetDuringTrialEstimate,4) < LowestDose |signif(TDtargetDuringTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",DuringTrialtargetprob*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetDuringTrialEstimate, colour="orange", lwd=1.1) +
annotate("text",label=paste(paste("TD",DuringTrialtargetprob*100),"Estimate"),
x=TDtargetDuringTrialEstimate,y=0,hjust=-0.1, vjust = -20,size=5,colour="orange")}
if (signif(TDtargetEndOfTrialEstimate,4) < LowestDose |signif(TDtargetEndOfTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",EndOfTrialtargetprob*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetEndOfTrialEstimate, colour="violet", lwd=1.1) +
annotate("text",label=paste(paste("TD",EndOfTrialtargetprob*100),"Estimate"),
x=TDtargetEndOfTrialEstimate,y=0,hjust=-0.1, vjust = -25,size=5,colour="violet")}
if (signif(Gstar,4) < LowestDose |signif(Gstar,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste("Estimated Gstar=",paste(Gstar," not within dose Grid")))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=Gstar, colour="green", lwd=1.1) +
annotate("text",label=" Gstar Estimate",
x=Gstar,y=0,hjust=-0.1, vjust = -25,size=5,colour="green")}
if (doselimit > max(data@doseGrid)){maxdoselimit<-max(data@doseGrid)} else {maxdoselimit <-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -35),
colour="red")
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0.1, vjust = -30),
colour="blue")
## return next best dose and plot
return(list(nextdose=ret,
DLEDuringTrialtarget=DuringTrialtargetprob,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
TDtargetDuringTrialAtDoseGrid=retD,
DLEEndOfTrialtarget=EndOfTrialtargetprob,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=retE,
GstarEstimate=Gstar,
GstarAtDoseGrid=Gstarret,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
CIGstar=CIGstar,
ratioGstar=ratioGstar,
plot=plot1))
} else if(is(Effmodel, "EffFlexi")) {
##first get the probDLE samples
points <- data@doseGrid
probDLESamples <- matrix(nrow=sampleSize(samples@options),
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:
probDLESamples[, i] <- prob(dose=points[i],
model=model,
samples=samples)
}
probDLE <- apply(probDLESamples,2,FUN=nextBest@TDderive)
DuringTrialtargetprob <- nextBest@DLEDuringTrialtarget
EndOfTrialtargetprob <- nextBest@DLEEndOfTrialtarget
## Find the TDtarget samples for During Trial and End of trial
TDtargetEndOfTrialSamples <- dose(prob=EndOfTrialtargetprob,
model=model,
samples=samples)
TDtargetDuringTrialSamples<-dose(prob=DuringTrialtargetprob,
model=model,
samples=samples)
## Find the TDtarget Estimate for During Trial and End of trial
TDtargetEndOfTrialEstimate <- as.numeric(nextBest@TDderive(TDtargetEndOfTrialSamples))
## Ensure the estimate is within dose range
#TDtargetEndOfTrialEstimate <- min(TDtargetEndOfTrialEstimate,max(data@doseGrid))
TDtargetDuringTrialEstimate<-as.numeric(nextBest@TDderive(TDtargetDuringTrialSamples))
## Ensure the estimate is within dose range
#TDtargetDuringTrialEstimate <- min(TDtargetDuringTrialEstimate,max(data@doseGrid))
##we have to get samples from the gain values at all dose levels
ExpEffsamples <- Effsamples@data$ExpEff
ExpEff <- apply(ExpEffsamples,2,FUN=nextBest@Gstarderive)
GainSamples <- matrix(nrow=sampleSize(samples@options),
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],
DLEmodel=model,
DLEsamples=samples,
Effmodel=Effmodel,
Effsamples=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]
##Obtain the Gstar estimate which is the 50th percentile of the Gstar samples
Gstar <- as.numeric(nextBest@Gstarderive(GstarSamples))
##Ensure the estimate is within dose range
#Gstar <- min(Gstar,max(data@doseGrid))
gainvalues <- apply(GainSamples,2,FUN=nextBest@Gstarderive)
dosesOK <- which(data@doseGrid <= doselimit)
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
##FIND the next dose which is the minimum between TDtargetDuringTrial and Gstar
nextdose<-min(TDtargetDuringTrialEstimate,Gstar)
##Find the dose level in doseGrid closest below nextdose
index <- suppressWarnings(max(which((signif(nextdose,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level in doseGrid closest below TDtargetEndOfTrial
indexE <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retE <- data@doseGrid[indexE]
##Find the dose level in doseGrid closest below TDtargetDuringTrial
indexD <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retD <- data@doseGrid[indexD]
##Find the dose level in doseGrid closest below Gstar
Gstarindex <- suppressWarnings(max(which((signif(Gstar,digits=4) - data@doseGrid[dosesOK]) >= 0)))
Gstarret <- data@doseGrid[Gstarindex]
##Find the 95% credibility interval of Gstar and its ratio of the upper to the lower limit
CIGstar <- as.numeric(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 <- as.numeric(quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975)))
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData<-data.frame(dose=rep(data@doseGrid,3),
values=c(probDLE,
ExpEff,
gainvalues))
gdata<-with(plotData,
data.frame(x=dose,
y=values,
group=c(rep("p(DLE)",length(data@doseGrid)),
rep("Expected Efficacy",length(data@doseGrid)),
rep("Gain",length(data@doseGrid))),
Type=factor("Estimate",levels="Estimate")
))
## produce plot
plot1 <- ggplot() +
geom_histogram(data=
data.frame(x=GstarSamples),
aes(x=x),
fill = "darkgreen", colour = "green3", binwidth = 25) +
xlab("Gstar")+ xlim(c(0,max(data@doseGrid)))+
ylab("Posterior density")
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
if (signif(TDtargetDuringTrialEstimate,4) < LowestDose |signif(TDtargetDuringTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",DuringTrialtargetprob*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plo1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetDuringTrialEstimate, colour="orange", lwd=1.1) +
annotate("text",label=paste(paste("TD",DuringTrialtargetprob*100),"Estimate"),
x=TDtargetDuringTrialEstimate,y=0,hjust=-0.1, vjust = -20,size=5,colour="orange")}
if (signif(TDtargetEndOfTrialEstimate,4) < LowestDose |signif(TDtargetEndOfTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",EndOfTrialtargetprob*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetEndOfTrialEstimate, colour="violet", lwd=1.1) +
annotate("text",label=paste(paste("TD",EndOfTrialtargetprob*100),"Estimate"),
x=TDtargetEndOfTrialEstimate,y=0,hjust=-0.1, vjust = -25,size=5,colour="violet")}
if (signif(Gstar,4) < LowestDose |signif(Gstar,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste("Estimated Gstar=",paste(Gstar," not within dose Grid")))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=Gstar, colour="green", lwd=1.1) +
annotate("text",label=" Gstar Estimate",
x=Gstar,y=0,hjust=+0.6, vjust = -25,size=5,colour="green")}
if (doselimit > max(data@doseGrid)){maxdoselimit<-max(data@doseGrid)} else {maxdoselimit <-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -35),
colour="red")
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0.1, vjust = -30),
colour="blue")
## return next best dose and plot
return(list(nextdose=ret,
DLEDuringTrialtarget=DuringTrialtargetprob,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
TDtargetDuringTrialAtDoseGrid=retD,
DLEEndOfTrialtarget=EndOfTrialtargetprob,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=retE,
GstarEstimate=Gstar,
GstarAtDoseGrid=Gstarret,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
CIGstar=CIGstar,
ratioGstar=ratioGstar,
plot=plot1))
} else stop("Effmodel needs to be of class Effloglog or EffFlexi")
})
## ------------------------------------------------------------------------------------------------
## 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("stopTrial",
signature=
signature(stopping="StoppingTDCIRatio",
dose="ANY",
samples="Samples",
model="ModelTox",
data="ANY"),
def=
function(stopping,dose,samples,model,data,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##check id targetEndOfTrial is a probability
stopifnot(is.probability(targetEndOfTrial))
## find the TDtarget End of Trial samples
TDtargetEndOfTrialSamples <- dose(prob=targetEndOfTrial,
model=model,
samples=samples)
##Find the upper and lower limit of the 95% credibility interval
CI <- quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975))
##The ratio of the upper to the lower 95% credibility interval
ratio <- as.numeric(CI[2]/CI[1])
##so can we stop?
doStop <- ratio <= stopping@targetRatio
##generate messgae
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"),
"targetRatio =", stopping@targetRatio)
##return both
return(structure(doStop,
messgae=text))
})
## ----------------------------------------------------------------------------------------------
## 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,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##check if targetEndOfTrial is a probability
stopifnot(is.probability(targetEndOfTrial))
## find the TDtarget End of Trial
TDtargetEndOfTrial <- dose(prob=targetEndOfTrial,
model=model)
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(targetEndOfTrial/(1-targetEndOfTrial))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval
CI <- c()
CI[2] <- exp(log(TDtargetEndOfTrial) + 1.96* sqrt(varEta))
CI[1] <- exp(log(TDtargetEndOfTrial) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratio <- as.numeric(CI[2]/CI[1])
##so can we stop?
doStop <- ratio <= stopping@targetRatio
##generate messgae
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"),
"targetRatio =", stopping@targetRatio)
##return both
return(structure(doStop,
messgae=text))
})
## --------------------------------------------------------------------------------------------------
## 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="StoppingGstarCIRatio",
dose="ANY",
samples="Samples",
model="ModelTox",
data="DataDual"),
def=
function(stopping,dose,samples,model,data,TDderive, Effmodel,Effsamples,Gstarderive,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##checks
stopifnot(is.probability(targetEndOfTrial))
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(prob=targetEndOfTrial,
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=sampleSize(samples@options),
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],
DLEmodel=model,
DLEsamples=samples,
Effmodel=Effmodel,
Effsamples=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@targetRatio
##generate messgae
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,"TDatrgetEndOfTrial estimate","Gstar estimate"), "is smaller with ratio =",
round(ratio,4), " which is " , ifelse(doStop,"is less than or equal to","greater than"),
"targetRatio =", stopping@targetRatio)
text <- c(text1,text2,text3)
##return both
return(structure(doStop,
message=text))
})
## -----------------------------------------------------------------------------------------------
## 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="StoppingGstarCIRatio",
dose="ANY",
samples="missing",
model="ModelTox",
data="DataDual"),
def=
function(stopping,dose,model,data,Effmodel,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##checks
stopifnot(is.probability(targetEndOfTrial))
stopifnot(is(Effmodel,"ModelEff"))
## find the TDtarget End of Trial
TDtargetEndOfTrial <- dose(prob=targetEndOfTrial,
model=model)
##Find the dose with maximum gain value
Gainfun<-function(DOSE){
-gain(DOSE,DLEmodel=model,Effmodel=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@c)} 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 <- t(deltaG)%*%covBETA%*%deltaG
##Find the upper and lower limit of the 95% credibility interval of Gstar
CIGstar <-c()
CIGstar[2] <- exp(logGstar + 1.96* sqrt(varlogGstar))
CIGstar[1] <- exp(logGstar - 1.96* sqrt(varlogGstar))
##The ratio of the upper to the lower 95% credibility interval
ratioGstar <- as.numeric(CIGstar[2]/CIGstar[1])
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(targetEndOfTrial/(1-targetEndOfTrial))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval of
##TDtargetEndOfTrial
CITDEOT <- c()
CITDEOT[2] <- exp(log(TDtargetEndOfTrial) + 1.96* sqrt(varEta))
CITDEOT[1] <- exp(log(TDtargetEndOfTrial) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(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@targetRatio
##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"),
"targetRatio =", stopping@targetRatio)
text <- c(text1,text2,text3)
##return both
return(structure(doStop,
message=text))
})
0liver0815/onc-crmpack-test documentation built on Feb. 19, 2022, 12:25 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.
#####################################################################################
## Author: Daniel Sabanes Bove [sabanesd *a*t* roche *.* com],
## Wai Yin Yeung [ w *.* yeung1 *a*t* lancaster *.* ac *.* uk]
## Project: Object-oriented implementation of CRM designs
##
## Time-stamp: <[Rules-methods.R] by DSB Die 09/06/2015 21:29>
##
## Description:
## Encapsulate the rule functions in formal methods.
##
## History:
## 07/02/2014 file creation
## 10/07/2015 Adding more classes
###################################################################################
##' @include Model-methods.R
##' @include Samples-class.R
##' @include Rules-class.R
{}
## ============================================================
## --------------------------------------------------
## Find out what is the next best dose
## --------------------------------------------------
##' Find the next best dose
##'
##' Compute the recommended next best dose.
##'
##' This function outputs the next best dose recommendation based on the
##' corresponding rule \code{nextBest}, the posterior \code{samples} from the
##' \code{model} and the underlying \code{data}.
##'
##' @param nextBest The rule, an object of class \code{\linkS4class{NextBest}}
##' @param doselimit The maximum allowed next dose. If this is an empty (length
##' 0) vector, then no dose limit will be applied in the course of dose
##' recommendation calculation, and a corresponding warning is given.
##' @param samples the \code{\linkS4class{Samples}} object
##' @param model The model input, an object of class \code{\linkS4class{Model}}
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots possible additional arguments without method dispatch
##' @return a list with the next best dose (element \code{value})
##' on the grid defined in \code{data}, and a plot depicting this recommendation
##' (element \code{plot}). In case of multiple plots also an element \code{singlePlots}
##' is included which returns the list of single plots, which allows for further
##' customization of these. Also additional list elements describing the outcome
##' of the rule can be contained.
##'
##' @export
##' @keywords methods
setGeneric("nextBest",
def=
function(nextBest, doselimit, samples, model, data, ...){
## there should be no default method,
## therefore just forward to next method!
standardGeneric("nextBest")
},
valueClass="list")
## --------------------------------------------------
## The MTD method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the MTD rule
##'
##' @example examples/Rules-method-NextBestMTD.R
##' @importFrom ggplot2 ggplot geom_density xlab ylab xlim aes geom_vline
##' geom_text
setMethod("nextBest",
signature=
signature(nextBest="NextBestMTD",
doselimit="numeric",
samples="Samples",
model="Model",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## First, generate the MTD samples.
mtdSamples <- dose(prob=nextBest@target,
model,
samples)
## then derive the next best dose
mtdEstimate <- nextBest@derive(mtdSamples=mtdSamples)
## be sure which doses are ok with respect to maximum
## possible dose - if one was specified
dosesOK <-
if(length(doselimit))
which(data@doseGrid <= doselimit)
else
seq_along(data@doseGrid)
## but now, round to the next possible grid point
index <- which.min(abs(data@doseGrid[dosesOK] - mtdEstimate))
ret <- data@doseGrid[dosesOK][index]
## produce plot
plot1 <- ggplot() +
geom_density(data=
data.frame(x=mtdSamples),
aes(x=x),
fill = "grey50", colour = "grey50") +
xlab("MTD") + ylab("Posterior density") +
xlim(range(data@doseGrid))
plot1 <- plot1 +
geom_vline(xintercept=mtdEstimate, colour="black", lwd=1.1) +
geom_text(data=
data.frame(x=mtdEstimate),
aes(x, 0,
label = "Est", hjust=+1, vjust = +1),
colour="black")
if(length(doselimit))
{
plot1 <- plot1 +
geom_vline(xintercept=doselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=doselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = +1),
colour="red")
}
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0, vjust = +1),
colour="blue")
## return next best dose and plot
return(list(value=ret,
plot=plot1))
})
## --------------------------------------------------
## The NCRM method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the NCRM method. The
##' additional list element \code{probs} contains the target and overdosing
##' probabilities (across all doses in the dose grid)
##' used in the derivation of the next best dose.
##'
##' @example examples/Rules-method-NextBestNCRM.R
##' @importFrom ggplot2 ggplot geom_bar xlab ylab ylim aes geom_vline
##' geom_hline geom_point
##' @importFrom gridExtra arrangeGrob
setMethod("nextBest",
signature=
signature(nextBest="NextBestNCRM",
doselimit="numeric",
samples="Samples",
model="Model",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## first we have to get samples from the dose-tox
## curve at the dose grid points.
probSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the probs, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
probSamples[, i] <- prob(dose=data@doseGrid[i],
model,
samples)
}
## Now compute probabilities to be in target and
## overdose tox interval
probTarget <-
colMeans((probSamples >= nextBest@target[1]) &
(probSamples <= nextBest@target[2]))
probOverdose <-
colMeans((probSamples > nextBest@overdose[1]) &
(probSamples <= nextBest@overdose[2]))
## which doses are eligible after accounting
## for maximum possible dose and discarding overdoses?
dosesBelowLimit <-
if(length(doselimit))
(data@doseGrid <= doselimit)
else
rep(TRUE, length(data@doseGrid))
dosesOK <- which(dosesBelowLimit &
(probOverdose < nextBest@maxOverdoseProb))
## check if there are doses that are OK
if(length(dosesOK))
{
## what is the recommended dose level?
## if maximum target probability is higher than some numerical
## threshold, then take that level, otherwise stick to the
## maximum level that is OK:
doseLevel <-
if(max(probTarget[dosesOK]) > 0.05)
{
which.max(probTarget[dosesOK])
} else {
which.max(data@doseGrid[dosesOK])
}
ret <- data@doseGrid[dosesOK][doseLevel]
} else {
## if none of the doses is OK:
doseLevel <- NA
ret <- NA
}
## produce plot
## first for the target probability
plot1 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probTarget * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="darkgreen",
fill="darkgreen") +
xlab("Dose") +
ylab(paste("Target probability [%]")) +
ylim(c(0, 100))
if(length(doselimit))
{
plot1 <- plot1 +
geom_vline(xintercept=doselimit,
lwd=1.1,
lty=2,
colour="black")
}
if(length(dosesOK))
{
plot1 <- plot1 +
geom_vline(xintercept=data@doseGrid[max(dosesOK)],
lwd=1.1,
lty=2,
colour="red")
plot1 <- plot1 +
geom_point(data=
data.frame(x=ret,
y=probTarget[dosesOK][doseLevel] *
100 + 0.03),
aes(x=x, y=y),
size=3,
pch=25,
col="red",
bg="red")
}
## second for the overdosing probability
plot2 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probOverdose * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="red",
fill="red") +
xlab("Dose") +
ylab("Overdose probability [%]") +
ylim(c(0, 100))
plot2 <- plot2 +
geom_hline(yintercept=nextBest@maxOverdoseProb * 100,
lwd=1.1,
lty=2,
colour="black")
## now plot them below each other
plotJoint <- gridExtra::arrangeGrob(plot1, plot2, nrow=2)
## return value and plot
return(list(value=ret,
plot=plotJoint,
singlePlots=list(plot1=plot1,
plot2=plot2),
probs=cbind(dose=data@doseGrid,
target=probTarget,
overdose=probOverdose)))
})
##' @describeIn nextBest Find the next best dose based on the NCRM method when
##' two parts trial is used - todo: need an example here for DataParts
##' @example examples/Rules-method-NextBestNCRM-DataParts.R
setMethod("nextBest",
signature=
signature(nextBest="NextBestNCRM",
doselimit="numeric",
samples="Samples",
model="Model",
data="DataParts"),
def=
function(nextBest, doselimit, samples, model, data, ...){
## exception when we are in part I or about to start part II!
if(all(data@part == 1L))
{
## here we will always propose the highest possible dose
## (assuming that the dose limit came from reasonable
## increments rule, i.e. inrementsRelativeParts)
if(identical(length(doselimit), 0L))
{
stop("doselimit needs to be given for Part I")
}
return(list(value=doselimit,
plot=NULL))
} else {
## otherwise we will just do the standard thing
callNextMethod(nextBest, doselimit, samples, model, data, ...)
}
})
## --------------------------------------------------
## The 3+3 method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the 3+3 method
##' @example examples/Rules-method-NextBestThreePlusThree.R
setMethod("nextBest",
signature=
signature(nextBest="NextBestThreePlusThree",
doselimit="missing",
samples="missing",
model="missing",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
## split the DLTs into the dose level groups
dltSplit <- split(data@y,
factor(data@x,
levels=data@doseGrid))
## number of patients and number of DLTs per dose level group
nPatients <- sapply(dltSplit, length)
nDLTs <- sapply(dltSplit, sum)
## what was the last dose level tested?
lastLevel <- tail(data@xLevel, 1)
## if there are less than 1/3 DLTs at that level
if(nDLTs[lastLevel]/nPatients[lastLevel] < 1/3)
{
## we could escalate, unless this is the highest
## level or the higher level was tried already
## (in which case it was not safe)
if((lastLevel == length(data@doseGrid)) ||
(nPatients[lastLevel+1] > 0))
{
nextLevel <- lastLevel
} else {
nextLevel <- lastLevel + 1
}
} else if(nDLTs[lastLevel]/nPatients[lastLevel] > 1/3) {
## rate here is too high, therefore deescalate
nextLevel <- lastLevel - 1
} else {
## otherwise: rate is 1/3 == 2/6,
## then it depends on the number of patients:
## if more than 3, then deescalate, otherwise stay.
nextLevel <-
if(nPatients[lastLevel] > 3)
lastLevel - 1
else
lastLevel
}
## do we stop here? only if we have no MTD
## or the next level has been tried enough (more than
## three patients already)
stopHere <-
if(nextLevel == 0)
{
TRUE
} else {
nPatients[nextLevel] > 3
}
## return value and plot
return(list(value=
if(nextLevel == 0) NA else data@doseGrid[nextLevel],
stopHere=stopHere))
})
## --------------------------------------------------
## The method for the dual endpoint model
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the dual endpoint
##' model. The additional list element \code{probs} contains the target and
##' overdosing probabilities (across all doses in the dose grid) used in the
##' derivation of the next best dose.
##' @example examples/Rules-method-NextBestDualEndpoint.R
##' @importFrom ggplot2 ggplot geom_bar xlab ylab ylim aes geom_vline
##' geom_hline geom_point
##' @importFrom gridExtra arrangeGrob
setMethod("nextBest",
signature=
signature(nextBest="NextBestDualEndpoint",
doselimit="numeric",
samples="Samples",
model="DualEndpoint",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## get the biomarker level samples
## at the dose grid points.
biomLevelSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the biomLevels, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
biomLevelSamples[, i] <- biomLevel(dose=data@doseGrid[i],
xLevel=i,
model,
samples)
}
## biomLevelSamples <- samples@data$betaW
## now get samples from the dose-tox
## curve at the dose grid points.
probSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the probs, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
probSamples[, i] <- prob(dose=data@doseGrid[i],
model,
samples)
}
## if target is relative to maximum
if(nextBest@scale == "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@data)){
## For each sample, look which dose is maximizing the
## simultaneous probability to be in the target biomarker
## range and below overdose toxicity
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1,ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= nextBest@target[1] * samples@data$Emax &
biomLevelSamples[, x] <= nextBest@target[2] * samples@data$Emax) /
nrow(biomLevelSamples)
})
}else{
## For each sample, look which was the minimum dose giving
## relative target level
targetIndex <- apply(biomLevelSamples, 1L,
function(x){
rnx <- range(x)
min(which((x >= nextBest@target[1] * diff(rnx) + rnx[1]) &
(x <= nextBest@target[2] * diff(rnx) + rnx[1] + 1e-10))
)
})
probTarget <- numeric(ncol(biomLevelSamples))
tab <- table(targetIndex)
probTarget[as.numeric(names(tab))] <- tab
probTarget <- probTarget / nrow(biomLevelSamples)
}
} else {
## otherwise 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
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1, ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= nextBest@target[1] &
biomLevelSamples[, x] <= nextBest@target[2]) /
nrow(biomLevelSamples)
})
}
## Now compute probabilities to be in
## overdose tox interval
probOverdose <-
colMeans((probSamples > nextBest@overdose[1]) &
(probSamples <= nextBest@overdose[2]))
## which doses are eligible after accounting
## for maximum possible dose and discarding overdoses?
dosesBelowLimit <-
if(length(doselimit))
(data@doseGrid <= doselimit)
else
rep(TRUE, length(data@doseGrid))
dosesOK <- which(dosesBelowLimit &
(probOverdose < nextBest@maxOverdoseProb))
## check if there are doses that are OK
if(length(dosesOK))
{
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
## what is the recommended dose level?
## if maximum target probability is higher than the numerical
## threshold, then take that level, otherwise stick to the
## maximum level that is OK:
doseLevel <-
if(max(probTarget[dosesOK]) > nextBest@targetThresh)
{
which.max(probTarget[dosesOK])
} else {
which.max(data@doseGrid[dosesOK])
}
ret <- data@doseGrid[dosesOK][doseLevel]
} else {
## if none of the doses is OK:
doseLevel <- NA
ret <- NA
}
## produce plot
## first for the target probability
plot1 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probTarget * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="darkgreen",
fill="darkgreen") +
xlab("Dose") +
ylab(paste("Target probability [%]")) +
ylim(c(0, 100))
if(length(doselimit))
{
plot1 <- plot1 +
geom_vline(xintercept=doselimit,
lwd=1.1,
lty=2,
colour="black")
}
if(length(dosesOK))
{
plot1 <- plot1 +
geom_vline(xintercept=data@doseGrid[max(dosesOK)],
lwd=1.1,
lty=2,
colour="red")
plot1 <- plot1 +
geom_point(data=
data.frame(x=ret,
y=probTarget[dosesOK][doseLevel] *
100 + 0.03),
aes(x=x, y=y),
size=3,
pch=25,
col="red",
bg="red")
}
## second for the overdosing probability
plot2 <- ggplot() +
geom_bar(data=
data.frame(x=data@doseGrid,
y=probOverdose * 100),
aes(x=x, y=y),
stat="identity",
position="identity",
width=min(diff(data@doseGrid)) / 2,
colour="red",
fill="red") +
xlab("Dose") +
ylab("Overdose probability [%]") +
ylim(c(0, 100))
plot2 <- plot2 +
geom_hline(yintercept=nextBest@maxOverdoseProb * 100,
lwd=1.1,
lty=2,
colour="black")
## now plot them below each other
plotJoint <- gridExtra::arrangeGrob(plot1, plot2, nrow=2)
## return value and plot
return(list(value=ret,
plot=plotJoint,
singlePlots=list(plot1=plot1,
plot2=plot2),
probs=cbind(dose=data@doseGrid,
target=probTarget,
overdose=probOverdose)))
})
## --------------------------------------------------
## The CRM MTD method
## --------------------------------------------------
##' @describeIn nextBest Find the next best dose based on CRM MTD allocation
##' probability
##'
##' @importFrom ggplot2 ggplot geom_density xlab ylab xlim aes geom_vline
##' geom_text
setMethod("nextBest",
signature=
signature(nextBest="NextBestMTDCRM",
doselimit="numeric",
samples="Samples",
model="Model",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
if(identical(length(doselimit), 0L))
{
warning("doselimit is empty, therefore no dose limit will be applied")
}
## first we have to get samples from the dose-tox
## curve at the dose grid points.
probSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the probs, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
probSamples[, i] <- prob(dose=data@doseGrid[i],
model,
samples)
}
## count the total number of samples
nrow <- sampleSize(samples@options)
## Now identify dose rate probabilities that are less than or
## equal to target
probTarget <-(probSamples <= nextBest@target)
## sum the number of doses for each row, that meet the
## criterion (TRUE)
sumMTD <- rowSums(probTarget[]=='TRUE')
## count the number of different values of the sums
countsumMTD <- table(sumMTD)
## if the dimension of countsumMTD ias not the same as the nGrid+1
## we need to impute with zero's the missing dose levels
if (dim(countsumMTD)<(data@nGrid+1)){
frame<-as.data.frame(countsumMTD)
framefull<-data.frame(sumMTD=0:data@nGrid)
frame1 <- merge(x=framefull ,y= frame,by="sumMTD", all.x = TRUE)
frame1$Freq[is.na(frame1$Freq)] = 0
countsumMTD <- as.table(frame1$Freq)
names(countsumMTD)<-c(0:data@nGrid)
}
## percentage of different values of the sums
freqMTD <- (countsumMTD/nrow)*100
## when the data object definition does NOT contain a placebo dose
if(data@placebo == "FALSE")
{
## if raw frequencies corresponding to planned doses are used for the nexbest dose
if(nextBest@pbomethod == "none")
{
countsumMTDnone<-countsumMTD[2:dim(countsumMTD)]
rr <- (countsumMTDnone/sum(countsumMTDnone))*100
rrmaxvalue<-rr[rr[]==max(rr)]
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## if cumulative frequencies corresponding to planned doses are used for the nexbest dose
## and the frequency of zero dose is added to the frequency of the lowest planned dose
if(nextBest@pbomethod == "min")
{
freqMTDdmin <-freqMTD
freqMTDdmin[2]<-freqMTDdmin[1]+freqMTDdmin[2]
rr<-freqMTDdmin[2:dim(freqMTDdmin)]
rrmaxvalue<-rr[rr[]==max(rr)]
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## if cumulative frequencies corresponding to planned doses are used for the nexbest dose
## and the frequency of zero dose is added to the frequency of the highest planned dose
if(nextBest@pbomethod == "max")
{
freqMTDdmax <-freqMTD
freqMTDdmax[dim(freqMTDdmax)]<-freqMTDdmax[1]+freqMTDdmax[dim(freqMTDdmax)]
rr<-freqMTDdmax[2:dim(freqMTDdmax)]
rrmaxvalue<-rr[rr[]==max(rr)]
}
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## when the data object definition contains a placebo dose
## if TRUE the first dose level in the grid is considered as PLACEBO
if(data@placebo == "TRUE")
{
rr<-freqMTD
rrmaxvalue<-rr[rr[]==max(rr)]
if (length(rrmaxvalue)==1){
doselevelMTD<-min(as.numeric(rownames(data.frame(rrmaxvalue))))
doseMTD<-data@doseGrid[doselevelMTD]
}
## in case of ties in rrmaxvalue's the lowest index/dose is selected
if (length(rrmaxvalue)>1){
doselevelMTD<-min(as.numeric(rownames(rrmaxvalue)))
doseMTD<-data@doseGrid[doselevelMTD]
}
}
## be sure which doses are ok with respect to maximum
## possible dose - if one was specified
dosesOK <-
if(length(doselimit)){
which(data@doseGrid <= doselimit)
}else{
seq_along(data@doseGrid)
}
dosesOKmax<-max(data@doseGrid[dosesOK])
if (doseMTD >= dosesOKmax){
ret <- dosesOKmax
}
if (doseMTD < dosesOKmax) {
ret <- doseMTD
}
## produce plot
plot1 <- ggplot() +
geom_line(data=data.frame(x=data@doseGrid,
y=as.numeric(rr) ),
aes(x=x, y=y),
linetype = "dashed",
size=1,
colour="red"
) +
geom_point(data=data.frame(x=data@doseGrid,
y=as.numeric(rr)
),
aes(x=x, y=y),
size=2,
colour="red"
) +
geom_text(data=data.frame(x=data@doseGrid,
y=as.numeric(rr)
),
mapping=aes(x=data@doseGrid,
y=as.numeric(rr) ,
label=paste(round(as.numeric(rr), 1),"%")
),
size = 3,
hjust = 0.5,
vjust = -1
)+
geom_vline(xintercept = doseMTD,
linetype = "dashed",
size=1,
colour="blue"
) +
annotate(geom = "text",
x = (doseMTD),
y = 97.5,
label = "Max",
colour="blue"
) +
geom_vline(xintercept = ret,
linetype = "dashed",
size=1,
colour="green") +
annotate(geom = "text",
x = ret,
y = 92.5,
colour="green",
label = "Next"
) +
xlab("Dose") +
ylab(paste("Allocation criterion [%]")) +
scale_x_continuous(limits=range(data@doseGrid) , breaks=data@doseGrid) +
scale_y_continuous(limits=c(0, 100) , breaks=seq(0,100,10))
## if doselimit exists add the dose range on the plot
if(length(doselimit))
{
plot1 <- plot1+
geom_rect(aes(xmin = data@doseGrid[1],
ymin = 0,
xmax =doselimit,
ymax = 100
),
alpha = 0.5
) +
annotate(geom = "text",
x = doselimit,
y = 95,
label = "Doselimit"
)
}
## return next best dose
return(list(value=ret,
plot=plot1))
})
## ============================================================
## --------------------------------------------------
## Determine the maximum possible next dose
## --------------------------------------------------
##' Determine the maximum possible next dose
##'
##' Determine the upper limit of the next dose based on the increments rule.
##'
##' This function outputs the maximum possible next dose, based on the
##' corresponding rule \code{increments} and the \code{data}.
##'
##' @param increments The rule, an object of class
##' \code{\linkS4class{Increments}}
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots further arguments
##' @return the maximum possible next dose
##'
##' @export
##' @keywords methods
setGeneric("maxDose",
def=
function(increments, data, ...){
## there should be no default method,
## therefore just forward to next method!
standardGeneric("maxDose")
},
valueClass="numeric")
## --------------------------------------------------
## The maximum allowable relative increments in intervals method
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' relative increments
##'
##' @example examples/Rules-method-maxDose-IncrementsRelative.R
setMethod("maxDose",
signature=
signature(increments="IncrementsRelative",
data="Data"),
def=
function(increments, data, ...){
## determine what was the last dose
lastDose <- tail(data@x, 1)
## determine in which interval this dose was
lastInterval <-
findInterval(x=lastDose,
vec=increments@intervals)
## so the maximum next dose is
ret <-
(1 + increments@increments[lastInterval]) *
lastDose
return(ret)
})
## --------------------------------------------------
## The maximum allowable number of dose levels method
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' maximum dose levels to increment for the next dose
##'
##' @example examples/Rules-method-maxDose-IncrementsNumDoseLevels.R
setMethod("maxDose",
signature=
signature(increments="IncrementsNumDoseLevels",
data="Data"),
def=
function(increments, data, ...){
## determine what was the level of the last dose
lastDoseLevel <- tail(data@xLevel, 1)
## determine the maximum next dose level
maxNextDoseLevel <- min(length(data@doseGrid),
lastDoseLevel + increments@maxLevels)
## so the maximum next dose is
ret <- data@doseGrid[maxNextDoseLevel]
return(ret)
})
## --------------------------------------------------
## The maximum allowable relative increments, with special rules for
## part 1 and beginning of part 2, method method
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' relative increments and part 1 and 2
##' @example examples/Rules-method-maxDose-IncrementsRelativeParts.R
setMethod("maxDose",
signature=
signature(increments="IncrementsRelativeParts",
data="DataParts"),
def=
function(increments, data, ...){
## determine if there are already cohorts
## belonging to part 2:
alreadyInPart2 <- any(data@part == 2L)
## if so, we just call the next higher method
if(alreadyInPart2)
{
callNextMethod(increments, data, ...)
} else {
## otherwise we have special rules.
## what dose level (index) has the highest dose
## so far?
lastDoseLevel <- matchTolerance(max(data@x),
data@part1Ladder)
## determine the next maximum dose
ret <-
if(data@nextPart == 1L)
{
## here the next cohort will still be in part 1.
## Therefore we just make one step on the part 1 ladder:
data@part1Ladder[lastDoseLevel + 1L]
} else {
## the next cohort will start part 2.
## if there was a DLT so far:
if(any(data@y == 1L))
{
data@part1Ladder[lastDoseLevel + increments@dltStart]
} else {
## otherwise
if(increments@cleanStart > 0)
{
## if we want to start part 2 higher than
## the last part 1 dose, use usual increments
callNextMethod(increments, data, ...)
} else {
## otherwise
data@part1Ladder[lastDoseLevel + increments@cleanStart]
}
}
}
return(ret)
}
})
## --------------------------------------------------
## The maximum allowable relative increments in terms of DLTs
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' relative increments determined by DLTs so far
##'
##' @example examples/Rules-method-maxDose-IncrementsRelativeDLT.R
setMethod("maxDose",
signature=
signature(increments="IncrementsRelativeDLT",
data="Data"),
def=
function(increments, data, ...){
## determine what was the last dose
lastDose <- tail(data@x, 1)
## determine how many DLTs have occurred so far
dltHappened <- sum(data@y)
## determine in which interval this is
interval <-
findInterval(x=dltHappened,
vec=increments@DLTintervals)
## so the maximum next dose is
ret <-
(1 + increments@increments[interval]) *
lastDose
return(ret)
})
## --------------------------------------------------
## The maximum allowable relative increments in terms of DLTs
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' multiple increment rules (taking the minimum across individual increments).
##'
##' @example examples/Rules-method-maxDose-IncrementMin.R
setMethod("maxDose",
signature=
signature(increments="IncrementMin",
data="Data"),
def=
function(increments, data, ...){
## apply the multiple increment rules
individualResults <-
sapply(increments@IncrementsList,
maxDose,
data=data,
...)
## so the maximum increment is the minimum across the individual increments
ret <- min(individualResults)
return(ret)
})
## --------------------------------------------------
## The maximum allowable number of dose levels method based on Hard Safety Rule - Beta(a,b)
## --------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' maximum eligible dose level to increment for the next dose
##' below the target toxicity, by a % probability based on the Beta(a,b) distribution.
##' this can be calculated for any number of patients in a cohort with a defined probability threshold of 90%
##' i.e. 1- pbeta(target, x+a, n-x+b) < 0.90
##' where x=DLTs, n=patients for a dose and use Beta(a=1,b=1) (or probably a=b=0.5)
##' if a dose satisfies the above criterio then eligible doses, are only the ones that that lower than that dose
##'
setMethod("maxDose",
signature=
signature(increments="IncrementsHSRBeta",
data="Data"),
def=
function(increments, data, ...){
#Load data dose vector and find the number of total subjects dosed at each dose
n<-table(data@x)
#Load data events vector at each dose and find the number of DLTs
yd<-table(data@y,data@x)
#check if both DLTs and no DLTs are present
if (nrow(yd)==2){
#Number of DLTs
y1<-yd[2,]
}
#check if only DLTs or only no DLTs are present
if (nrow(yd)==1){
#Number of DLTs in each scenario:
#if only DLTs
if (dimnames(yd)[1]==1){y1=n}
#if only no DLTs
if (dimnames(yd)[1]==0){y1=0}
}
#Calculate the probability of DLTs at each dose tested
# i.e. 1- pbeta(target, y+a, n-y+b) < 0.90
# where y=DLTs, n=patients for a dose and use Beta(a,b)
probTarget<-(1- pbeta(increments@target, y1+increments@a, n-y1+increments@b))
# find toxic doses
toxic<-probTarget>=increments@prob
# find safe doses
safe<-probTarget<increments@prob
# are there any safe doses
safegrid<-sum(safe[]=="TRUE")
# are there any toxic doses
toxicgrid<-sum(toxic[]=="TRUE")
# if only safe doses exist (toxic doses do not exist all) then all grid doses are eligible and max next dose is the maximum planned dose grid
if (safegrid[]>0 && toxicgrid[]==0){
maxdose<-max(data@doseGrid[])
}else if (safegrid[]>0 && toxicgrid[]>0){
# if safe and toxic doses both exist, then max next dose is the maximum safe dose
maxdose<-data@doseGrid[safegrid]
}else if (safegrid[]==0 && toxicgrid[]>0){
# if safe doses do not exist and all doses toxic
maxdose<-min(data@doseGrid)
# no dose is eligible so we assign the lowest (first) dose. This will force the next recommended dose to be the lowest dose and will trigger the stopping rule stoppingruleHSFBeta.
# therefore when using this increment rule it is necessary that you also use the relevant stoppingruleHSFBeta with the same setting. In this
# way the absence of safe doses including the lowest one will trigger trial to stop as long as these 2 are set the same.
}
## so the maximum next dose is
ret <- maxdose
return(ret)
})
## -----------------------------------------------------------------------------------------
## The maximum allowable number of dose levels method relative to the max dose tested so far
## -----------------------------------------------------------------------------------------
##' @describeIn maxDose Determine the maximum possible next dose based on
##' maximum dose levels to increment for the next dose
##'
setMethod("maxDose",
signature=
signature(increments="IncrementsNumDoseLevelsMaxTested",
data="Data"),
def=
function(increments, data, ...){
## determine what was the level of the max dose tested so far
maxDoseLevel <- max(data@xLevel)
## determine the maximum next dose level
maxNextDoseLevel <- min(length(data@doseGrid),
maxDoseLevel + increments@maxLevels)
## so the maximum next dose is
ret <- data@doseGrid[maxNextDoseLevel]
return(ret)
})
## ============================================================
## --------------------------------------------------
## "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@stopList <- c(e1@stopList,
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@stopList <- c(e1,
e2@stopList)
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@stopList <- c(e1@stopList,
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@stopList <- c(e1,
e2@stopList)
return(e2)
})
## --------------------------------------------------
## Stop the trial?
## --------------------------------------------------
##' Stop the trial?
##'
##' This function returns whether to stop the trial.
##'
##' @param stopping The rule, an object of class
##' \code{\linkS4class{Stopping}}
##' @param dose the recommended next best dose
##' @param samples the \code{\linkS4class{Samples}} object
##' @param model The model input, an object of class \code{\linkS4class{Model}}
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots additional arguments
##'
##' @return logical value: \code{TRUE} if the trial can be stopped, \code{FALSE}
##' otherwise. It should have an attribute \code{message} which gives the reason
##' for the decision.
##'
##' @export
##' @example examples/Rules-method-CombiningStoppingRulesAndOr.R
##' @keywords methods
setGeneric("stopTrial",
def=
function(stopping, dose, samples, model, data, ...){
## if the recommended next dose is NA,
## stop in any case.
if(is.na(dose))
{
return(structure(TRUE,
message="Recommended next best dose is NA"))
} else if(data@placebo && dose == min(data@doseGrid))
{
return(structure(TRUE,
message="Recommended next best dose is placebo dose"))
}
## there should be no default method,
## therefore just forward to next method!
standardGeneric("stopTrial")
},
valueClass="logical")
## --------------------------------------------------
## 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@stopList,
stopTrial,
dose=dose,
model=model,
data=data,
...)
} else {
lapply(stopping@stopList,
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))
})
## --------------------------------------------------
## 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@stopList,
stopTrial,
dose=dose,
model=model,
data=data,
...)
} else {
lapply(stopping@stopList,
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))
})
## --------------------------------------------------
## 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@stopList,
stopTrial,
dose=dose,
model=model,
data=data,
...)
} else {
lapply(stopping@stopList,
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))
})
## --------------------------------------------------
## 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))
})
## -------------------------------------------------------------
## 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 <- 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))
})
## --------------------------------------------------
## 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))
})
## --------------------------------------------------
## 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))
})
## --------------------------------------------------
## Stopping based on probability of target tox interval
## --------------------------------------------------
##' @describeIn stopTrial Stop based on probability of target tox interval
##'
##' @example examples/Rules-method-stopTrial-StoppingTargetProb.R
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetProb",
dose="numeric",
samples="Samples",
model="Model",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## first we have to get samples from the dose-tox
## curve at the dose.
probSamples <- prob(dose=dose,
model,
samples)
## Now compute probability to be in target interval
probTarget <-
mean((probSamples >= stopping@target[1]) &
(probSamples <= stopping@target[2]))
## so can we stop?
doStop <- probTarget >= stopping@prob
## generate message
text <-
paste("Probability for target toxicity is",
round(probTarget * 100),
"% for dose",
dose,
"and thus",
ifelse(doStop, "above", "below"),
"the required",
round(stopping@prob * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## --------------------------------------------------
## 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="Model",
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(prob=stopping@target,
model,
samples)
## what is the absolute threshold?
absThresh <- stopping@thresh * dose
## what is the probability to be above this dose?
prob <- 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, "above", "below"),
"the required",
round(stopping@prob * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## --------------------------------------------------
## Stopping based on probability of targeting biomarker
## --------------------------------------------------
##' @describeIn stopTrial Stop based on probability of targeting biomarker
##'
##' @example examples/Rules-method-stopTrial-StoppingTargetBiomarker.R
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetBiomarker",
dose="numeric",
samples="Samples",
model="DualEndpoint",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## compute the target biomarker prob at this dose
## get the biomarker level samples
## at the dose grid points.
biomLevelSamples <- matrix(nrow=sampleSize(samples@options),
ncol=data@nGrid)
## evaluate the biomLevels, for all samples.
for(i in seq_len(data@nGrid))
{
## Now we want to evaluate for the
## following dose:
biomLevelSamples[, i] <- biomLevel(dose=data@doseGrid[i],
xLevel=i,
model,
samples)
}
## if target is relative to maximum
if(stopping@scale == "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@data)){
## For each sample, look which dose is maximizing the
## simultaneous probability to be in the target biomarker
## range and below overdose toxicity
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1,ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= stopping@target[1]*samples@data$Emax &
biomLevelSamples[, x] <= stopping@target[2]*samples@data$Emax) / nrow(biomLevelSamples)
})
}else{
## For each sample, look which was the minimum dose giving
## relative target level
targetIndex <- apply(biomLevelSamples, 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))
)
})
probTarget <- numeric(ncol(biomLevelSamples))
tab <- table(targetIndex)
probTarget[as.numeric(names(tab))] <- tab
probTarget <- probTarget / nrow(biomLevelSamples)
}
} else {
## otherwise 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
probTarget <- numeric(ncol(biomLevelSamples))
probTarget <- sapply(seq(1, ncol(biomLevelSamples)),
function(x){
sum(biomLevelSamples[, x] >= stopping@target[1] &
biomLevelSamples[, x] <= stopping@target[2]) /
nrow(biomLevelSamples)
})
}
## so for this dose we have:
probTarget <- probTarget[which(data@doseGrid == dose)]
## so can we stop?
doStop <- probTarget >= stopping@prob
## generate message
text <-
paste("Probability for target biomarker is",
round(probTarget * 100),
"% for dose",
dose,
"and thus",
ifelse(doStop, "above", "below"),
"the required",
round(stopping@prob * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## --------------------------------------------------
## 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 <- (dose == data@doseGrid[data@nGrid])
return(structure(isHighestDose,
message=
paste("Next best dose is", dose, "and thus",
ifelse(isHighestDose, "the",
"not the"),
"highest dose")))
})
## -------------------------------------------------------------
## Stopping based on the lowest dose meeting the hard safety criteria using Beta based DLT probability
## -------------------------------------------------------------
##' @describeIn stopTrial Stop based on the lowest dose meeting the hard
##' safety criteria using Beta based DLT probability
setMethod("stopTrial",
signature=
signature(stopping="StoppingLowestDoseHSRBeta",
dose="numeric",
samples="ANY",
model="ANY",
data="Data"),
def=
function(stopping, dose, samples, model, data, ...){
#Load data dose vector and find the number of total subjects dosed at this the lowest planned dose
dose <- data@doseGrid[1]
if (dose %in% data@x){
xd<-table(data@x==dose)
if (dim(xd)>0){
n<-xd["TRUE"]
#Load data events vector at this dose and find the number of DLTs
yd<-table(data@y,data@x==dose)
#check if both DLTs and no DLTs are present
if (nrow(yd)==2){
#Number of DLTs
y1<-yd[2,"TRUE"]
}
#check if only DLTs or only no DLTs are present
if (nrow(yd)==1){
#Number of DLTs in each scenario:
#if only DLTs
if (dimnames(yd)[1]==1){y1=xd}
#if only no DLTs
if (dimnames(yd)[1]==0){y1=0}
}
#Calculate the probability of DLTs
# i.e. 1- pbeta(target, y+a, n-y+b) < 0.90
# where y=DLTs, n=patients for a dose and use Beta(a,b)
probTarget<-(1- pbeta(stopping@target, y1+stopping@a, n-y1+stopping@b))
## so can we stop?
doStop <- probTarget >= stopping@prob
## generate message
text <-
paste("Probability that the lowest dose",dose,"being toxic based on Beta(",
stopping@a,",",stopping@b,") is",probTarget,
"and thus",
ifelse(doStop, "above", "below"),
"the prespecified risk of",
stopping@prob)
## return both
return(structure(doStop,
message=text))
}
}else{
probTarget<-NA
## so can we stop?
doStop <- is.numeric(probTarget)
## generate message
text <-
paste("Lowest dose",dose,"never tested before and thus probability the lowest dose being toxic based on Beta(",
stopping@a,",",stopping@b,") is",probTarget
)
## return both
return(structure(doStop,
message=text))
}
})
## --------------------------------------------------
## Stopping based on precision of MTD calculated as CV(MTD)
## --------------------------------------------------
##' @describeIn stopTrial Stop based on precision of MTD calculated as CV(MTD)
setMethod("stopTrial",
signature=
signature(stopping="StoppingMTDCV",
dose="numeric",
samples="Samples",
model="Model",
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(prob=stopping@target,
model,
samples)
## CV of MTD derived based on MTD samples
MTDCV <- (mad(MTDSamples)/median(MTDSamples))
## so can we stop?
doStop <- ((MTDCV <= stopping@thresh) & (MTDCV >= 0))
## generate message
text <-
paste("CV of MTD is",
round(MTDCV * 100),"% and thus",
ifelse(doStop, "below", "above"),
"the required precision threshold of",
round(stopping@thresh * 100),
"%")
## return both
return(structure(doStop,
message=text))
})
## -----------------------------------------------------------------------------------------------
## Stopping based on probability of target tox interval and number of patients near to the highest dose
## -----------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on probability of target tox interval and
##' number of patients near to the highest dose
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetProbPatientsNearHighestDose",
dose="numeric",
samples="Samples",
model="Model",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## first we have to get samples from the dose-tox
## curve at the dose.
dose <- max(data@doseGrid)
probSamples <- prob(dose=dose,
model,
samples)
## Now compute probability to be in target interval
probTarget <-
mean((probSamples >= stopping@target[1]) &
(probSamples <= stopping@target[2]))
## so can we stop due to toxicity interval?
doStopprobTarget <- probTarget >= stopping@prob
lower <- (100 - stopping@percentage) / 100 * dose
upper <- (100 + stopping@percentage) / 100 * dose
## how many patients lie there?
nPatients <- sum((data@x >= lower) & (data@x <= upper))
## so can we stop due to sufficient number of patients?
doStopnPatients <- nPatients >= stopping@nPatients
## so can we stop overall?
doStop <- all(c(doStopprobTarget,doStopnPatients)=="TRUE")
## generate message
text <-
paste("Probability for target toxicity is ",
round(probTarget * 100),
" % for highest dose ",
dose,
" and thus ",
ifelse(doStopprobTarget, "above", "below"),
" the required ",
round(stopping@prob * 100),
"%"," and ",
nPatients,
" patients lie within ",
stopping@percentage,
"% of the highest dose ",
dose,
". This ",
ifelse(doStopnPatients, "reached", "is below"),
" the required ",
stopping@nPatients,
" patients",
sep="")
## return both
return(structure(doStop,
message=text))
})
## -----------------------------------------------------------------------------------------------
## Stopping based on probability of target tox interval and number of patients near to the lowest dose
## -----------------------------------------------------------------------------------------------
##' @describeIn stopTrial Stop based on probability of target tox interval and
##' number of patients near to the lowest dose
setMethod("stopTrial",
signature=
signature(stopping="StoppingTargetProbPatientsNearLowestDose",
dose="numeric",
samples="Samples",
model="Model",
data="ANY"),
def=
function(stopping, dose, samples, model, data, ...){
## first we have to get samples from the dose-tox
## curve at the dose.
dose <- data@doseGrid[1]
probSamples <- prob(dose=dose,
model,
samples)
## Now compute probability to be in target interval
probTarget <-
mean((probSamples >= stopping@target[1]) &
(probSamples <= stopping@target[2]))
## so can we stop due to toxicity interval?
doStopprobTarget <- probTarget >= stopping@prob
lower <- (100 - stopping@percentage) / 100 * dose
upper <- (100 + stopping@percentage) / 100 * dose
## how many patients lie there?
nPatients <- sum((data@x >= lower) & (data@x <= upper))
## so can we stop due to sufficient number of patients?
doStopnPatients <- nPatients >= stopping@nPatients
## so can we stop overall?
doStop <- all(c(doStopprobTarget,doStopnPatients)=="TRUE")
## generate message
text <-
paste("Probability for target toxicity is",
round(probTarget * 100),
" % for lowest dose ",
dose,
" and thus ",
ifelse(doStopprobTarget, "above", "below"),
" the required ",
round(stopping@prob * 100),
"%"," and ",
nPatients,
" patients lie within ",
stopping@percentage,
" % of the lowest dose ",
dose,
". This ",
ifelse(doStopnPatients, "reached", "is below"),
" the required ",
stopping@nPatients,
" patients",
sep="")
## return both
return(structure(doStop,
message=text))
})
## ============================================================
## --------------------------------------------------
## "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(...))
})
## --------------------------------------------------
## Determine the size of the next cohort
## --------------------------------------------------
##' Determine the size of the next cohort
##'
##' This function determines the size of the next cohort.
##'
##' @param cohortSize The rule, an object of class
##' \code{\linkS4class{CohortSize}}
##' @param dose the next dose
##' @param data The data input, an object of class \code{\linkS4class{Data}}
##' @param \dots additional arguments
##'
##' @return the size as integer value
##'
##' @export
##' @keywords methods
setGeneric("size",
def=
function(cohortSize, dose, data, ...){
## if the recommended next dose is NA,
## don't check and return 0
if(is.na(dose))
{
return(0L)
}
## there should be no default method,
## therefore just forward to next method!
standardGeneric("size")
},
valueClass="integer")
## --------------------------------------------------
## The dose range method
## --------------------------------------------------
##' @describeIn size Determine the cohort size based on the range into which the
##' next dose falls into
##'
##' @example examples/Rules-method-size-CohortSizeRange.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeRange",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## determine in which interval the next dose is
interval <-
findInterval(x=dose,
vec=cohortSize@intervals)
## so the cohort size is
ret <- cohortSize@cohortSize[interval]
return(ret)
})
## --------------------------------------------------
## The DLT range method
## --------------------------------------------------
##' @describeIn size Determine the cohort size based on the number of DLTs so
##' far
##'
##' @example examples/Rules-method-size-CohortSizeDLT.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeDLT",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## determine how many DLTs have occurred so far
dltHappened <- sum(data@y)
## determine in which interval this is
interval <-
findInterval(x=dltHappened,
vec=cohortSize@DLTintervals)
## so the cohort size is
ret <- cohortSize@cohortSize[interval]
return(ret)
})
## --------------------------------------------------
## Size based on maximum of multiple cohort size rules
## --------------------------------------------------
##' @describeIn size Size based on maximum of multiple cohort size rules
##' @example examples/Rules-method-size-CohortSizeMax.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeMax",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## evaluate the individual cohort size rules
## in the list
individualResults <-
sapply(cohortSize@cohortSizeList,
size,
dose=dose,
data=data,
...)
## summarize to obtain overall result
overallResult <- max(individualResults)
return(overallResult)
})
## --------------------------------------------------
## Size based on minimum of multiple cohort size rules
## --------------------------------------------------
##' @describeIn size Size based on minimum of multiple cohort size rules
##' @example examples/Rules-method-size-CohortSizeMin.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeMin",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
## evaluate the individual cohort size rules
## in the list
individualResults <-
sapply(cohortSize@cohortSizeList,
size,
dose=dose,
data=data,
...)
## summarize to obtain overall result
overallResult <- min(individualResults)
return(overallResult)
})
## --------------------------------------------------
## Constant cohort size
## --------------------------------------------------
##' @describeIn size Constant cohort size
##' @example examples/Rules-method-size-CohortSizeConst.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeConst",
dose="ANY",
data="Data"),
def=
function(cohortSize, dose, data, ...){
return(cohortSize@size)
})
## --------------------------------------------------
## Cohort size based on the parts
## --------------------------------------------------
##' @describeIn size Cohort size based on the parts
##' @example examples/Rules-method-size-CohortSizeParts.R
setMethod("size",
signature=
signature(cohortSize="CohortSizeParts",
dose="ANY",
data="DataParts"),
def=
function(cohortSize, dose, data, ...){
return(cohortSize@sizes[data@nextPart])
})
## ================================================================
## The nextBest method based only on DLE responses with samples
## ================================================================
##' @describeIn nextBest Find the next best dose based on the 'NextBestTDsamples'
##' class rule. This a method based only on the DLE responses and for
##' \code{\linkS4class{LogisticIndepBeta}} model class object involving DLE samples
##'
##' @importFrom ggplot2 ggplot geom_density xlab ylab xlim aes geom_vline
##' geom_text
##'
##' @example examples/Rules-method-nextbest_TDsamples.R
##'
##' @export
##' @keywords methods
setMethod("nextBest",
signature=
signature(nextBest="NextBestTDsamples",
doselimit="numeric",
samples="Samples",
model="LogisticIndepBeta",
data="Data"),
def=
function(nextBest, doselimit, samples, model, data, ...){
## First, generate the TDtargetDuringTrial (TDtarget During a Trial) and
## TDtargetEndOfTrial (TDtarget at the EndOfTrial) samples.
TDtargetDuringTrialSamples <- dose(prob=nextBest@targetDuringTrial,
model,
samples)
TDtargetEndOfTrialSamples <- dose(prob=nextBest@targetEndOfTrial,
model,
samples)
## then derive the prior/posterior mean of the above two samples
TDtargetDuringTrialEstimate <- as.numeric(nextBest@derive(TDsamples=TDtargetDuringTrialSamples))
TDtargetEndOfTrialEstimate <- as.numeric(nextBest@derive(TDsamples=TDtargetEndOfTrialSamples))
## be sure which doses are ok with respect to maximum
## possible dose
dosesOK <- which(data@doseGrid <= doselimit)
##Find the index of next dose in the doseGrid
##next dose is the dose level closest below the TDtargetDuringTrialEstimate
index <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level (in doseGrid) closest below the TDtargetEndOfTrialEstimate
index1 <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret1 <- data@doseGrid[dosesOK][index1]
CITDEOT <- as.numeric(quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975)))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
## produce plot
plot1 <- ggplot() +
geom_density(data=
data.frame(x=TDtargetDuringTrialSamples),
aes(x=x),
fill = "grey50", colour = "grey50")
plot1 <- plot1 +
geom_density(data=
data.frame(x=TDtargetEndOfTrialSamples),
aes(x=x),
fill = "grey50", colour = "violet") +
xlab("TD") + ylab("Posterior density") +
xlim(range(data@doseGrid))
plot1 <- plot1+
geom_vline(xintercept=TDtargetDuringTrialEstimate, colour="orange", lwd=1.1) +
annotate("text",label=paste(paste("TD",nextBest@targetDuringTrial*100),"Estimate"),
x=TDtargetDuringTrialEstimate,y=0,hjust=-0.1, vjust = -20,size=5,colour="orange")
plot1 <- plot1+
geom_vline(xintercept=TDtargetEndOfTrialEstimate, colour="violet", lwd=1.1) +
annotate("text",label=paste(paste("TD",nextBest@targetEndOfTrial*100),"Estimate"),
x=TDtargetEndOfTrialEstimate,y=0,hjust=-0.1, vjust = -25,size=5,colour="violet")
if (doselimit > max(data@doseGrid)){maxdoselimit<-max(data@doseGrid)} else {maxdoselimit <-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -35),
colour="red")
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0.1, vjust = -30),
colour="blue")
## return next best dose and plot
return(list(nextdose=ret,
targetDuringTrial=nextBest@targetDuringTrial,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
targetEndOfTrial=nextBest@targetEndOfTrial,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=ret1,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
plot=plot1))
})
## -------------------------------------------------------------------------------
## The nextBest method based only on DLE responses without samples
## -----------------------------------------------------------------------------
##' @describeIn nextBest Find the next best dose based on the 'NextBestTD'
##' class rule. This a method based only on the DLE responses and for
##' \code{\linkS4class{LogisticIndepBeta}} model class object without DLE samples
##'
##' @param SIM internal command to notify if this method is used within simulations. Default as FALSE
##' @importFrom ggplot2 ggplot xlab ylab xlim aes geom_vline
##' geom_text
##'
##' @example examples/Rules-method-nextbest_TD.R
##'
##' @export
##' @keywords methods
setMethod("nextBest",
signature=
signature(nextBest="NextBestTD",
doselimit="numeric",
samples="missing",
model="LogisticIndepBeta",
data="Data"),
def=
function(nextBest, doselimit, model, data, SIM=FALSE,...){
##Find the target prob During Trial
targetDuringTrial<- nextBest@targetDuringTrial
##Find the target prob End of Trial
targetEndOfTrial<- nextBest@targetEndOfTrial
mylabel<-targetDuringTrial*100
label2 <-targetEndOfTrial*100
## Find the TD30 Estimate and TD(target) Estimate
TDtargetEndOfTrialEstimate <- dose(prob=targetEndOfTrial,
model)
TDEfourdg<-signif(TDtargetEndOfTrialEstimate,digits=4)
TDtargetDuringTrialEstimate<-dose(prob=targetDuringTrial,
model)
TDDfourdg<-signif(TDtargetDuringTrialEstimate,digits=4)
probDLE=prob(dose=data@doseGrid,
model=model)
## be sure which doses are ok with respect to maximum
## possible dose
dosesOK <- which(data@doseGrid <= doselimit)
##Find the index of next dose in the doseGrid
##next dose is the dose level closest below the TDtargetEstimate
index <- suppressWarnings(max(which((TDDfourdg- data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level (in doseGrid) closest below the TD30Estimate
index <- suppressWarnings(max(which((TDEfourdg - data@doseGrid[dosesOK]) >= 0)))
retTDE <- data@doseGrid[dosesOK][index]
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(targetEndOfTrial/(1-targetEndOfTrial))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval
CITDEOT <- c()
CITDEOT[2] <- exp(log(TDtargetEndOfTrialEstimate) + 1.96* sqrt(varEta))
CITDEOT[1] <- exp(log(TDtargetEndOfTrialEstimate) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData <- data.frame(dose=data@doseGrid,
probDLE=prob(dose=data@doseGrid,
model=model))
##make the plot
gdata <- with(plotData,
data.frame(x=dose,
y=probDLE,
group=rep("Estimated DLE",each=nrow(plotData)),
Type=factor(rep("Estimated DLE",nrow(plotData)),levels="Estimated DLE")))
plot1 <- ggplot(data=gdata, aes(x=x,y=y), group=group) +
xlab("Dose Levels")+
ylab(paste("Probability of DLE")) + ylim(c(0,1)) + xlim(c(0,max(data@doseGrid))) +
geom_line(colour=I("red"), size=1.5)
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
if ((TDDfourdg < LowestDose)|(TDDfourdg > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("TD",targetDuringTrial*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plot1<- plot1}
} else {plot1 <- plot1+
geom_point(data=data.frame(x=TDtargetDuringTrialEstimate,y=targetDuringTrial),aes(x=x,y=y),colour="orange", shape=15, size=8) +
annotate("text",label=paste(paste("TD",mylabel),"Estimate"),x=TDtargetDuringTrialEstimate+1,y=targetDuringTrial-0.2,size=5,colour="orange")}
if ((TDEfourdg < LowestDose)|(TDEfourdg > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("TD",targetEndOfTrial*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {plot1 <- plot1+
geom_point(data=data.frame(x=TDtargetEndOfTrialEstimate,y=0.3),aes(x=x,y=y),colour="violet", shape=16, size=8) +
annotate("text",label=paste(paste("TD",label2),"Estimate"),x=TDtargetEndOfTrialEstimate+1,y=targetEndOfTrial-0.1,size=5,colour="violet")}
if (doselimit > max(data@doseGrid)) {maxdoselimit <- max(data@doseGrid)} else {maxdoselimit<-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="brown", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -30),
colour="brown")
plot1 <-plot1 +
geom_vline(xintercept=ret, colour="purple", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0, vjust = -30),
colour="purple")
## return next best dose and plot
return(list(nextdose=ret,
targetDuringTrial=targetDuringTrial,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
targetEndOfTrial=targetEndOfTrial,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialatdoseGrid=retTDE,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
plot=plot1))
})
## ------------------------------------------------------------------------------------
## the nextBest method based on DLE and efficacy responses without DLE and efficacy samples
## -------------------------------------------------------------------------- ----------
##' @describeIn nextBest for slots \code{nextBest},\code{doselimit}, \code{data} and \code{SIM}. This is
##' a function to find the next best dose based on the 'NextBestMaxGain'
##' class rule. This a method based on the DLE responses and efficacy responses without DLE and
##' efficacy samples.
##'
##' @param Effmodel the efficacy model of \code{\linkS4class{ModelEff}} class object
##'
##' @importFrom ggplot2 ggplot xlab ylab xlim aes geom_vline
##' geom_text
##'
##' @example examples/Rules-method-nextbest_MaxGain.R
##'
##' @importFrom ggplot2 scale_colour_manual
##' @export
##' @keywords methods
setMethod("nextBest",
signature=
signature(nextBest="NextBestMaxGain",
doselimit="numeric",
samples="missing",
model="ModelTox",
data="DataDual"),
def=
function(nextBest,doselimit,model,data,Effmodel,SIM=FALSE,...){
stopifnot(is(Effmodel, "ModelEff"))
DuringTrialtargetprob <- nextBest@DLEDuringTrialtarget
EndOfTrialtargetprob <- nextBest@DLEEndOfTrialtarget
## Find the TDtarget Estimate for During Trial and End of trial
TDtargetEndOfTrialEstimate <- dose(prob=EndOfTrialtargetprob,model=model)
TDtargetDuringTrialEstimate<-dose(prob=DuringTrialtargetprob,model=model)
##Get all prob of DLE at all dose levels
probDLE=prob(dose=data@doseGrid,
model=model)
##Define gain function
Gainfun<-function(DOSE){
-gain(DOSE,DLEmodel=model,Effmodel=Effmodel)
}
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
##Find the dose which gives the maximum gain
Gstar<-(optim(LowestDose,Gainfun, method = "L-BFGS-B", lower=LowestDose,upper=max(data@doseGrid))$par)
##Find the maximum gain value
MaxGain<--(optim(LowestDose,Gainfun,method = "L-BFGS-B", lower=LowestDose,upper=max(data@doseGrid))$value)
## be sure which doses are ok with respect to maximum
## possible dose
dosesOK <- which(data@doseGrid <= doselimit)
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
##FIND the next dose which is the minimum between TDtargetDuringTrial and Gstar
nextdose<-min(TDtargetDuringTrialEstimate,Gstar)
##Find the dose level in doseGrid closest below nextdose
index <- suppressWarnings(max(which((signif(nextdose,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level in doseGrid closest below TDtargetEndOfTrial
indexE <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retE <- data@doseGrid[indexE]
##Find the dose level in doseGrid closest below TDtargetDuringTrial
indexD <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retD <- data@doseGrid[indexD]
##Find the dose level in doseGrid closest below Gstar
Gstarindex <- suppressWarnings(max(which((signif(Gstar,digits=4) - data@doseGrid[dosesOK]) >= 0)))
Gstarret <- data@doseGrid[Gstarindex]
if (data@placebo){
logGstar <- log(Gstar+Effmodel@c)} 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
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 <- t(deltaG)%*%covBETA%*%deltaG
##Find the upper and lower limit of the 95% credibility interval of Gstar
CIGstar <-c()
CIGstar[2] <- exp(logGstar + 1.96* sqrt(varlogGstar))
CIGstar[1] <- exp(logGstar - 1.96* sqrt(varlogGstar))
##The ratio of the upper to the lower 95% credibility interval
ratioGstar <- as.numeric(CIGstar[2]/CIGstar[1])
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(EndOfTrialtargetprob/(1-EndOfTrialtargetprob))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval of
##TDtargetEndOfTrial
CITDEOT <- c()
CITDEOT[2] <- exp(log(TDtargetEndOfTrialEstimate) + 1.96* sqrt(varEta))
CITDEOT[1] <- exp(log(TDtargetEndOfTrialEstimate) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData<-data.frame(dose=rep(data@doseGrid,3),
values=c(prob(dose=data@doseGrid,
model=model),
ExpEff(dose=data@doseGrid,
model=Effmodel),
gain(dose=data@doseGrid,
DLEmodel=model,
Effmodel=Effmodel)))
gdata<-with(plotData,
data.frame(x=dose,
y=values,
group=c(rep("p(DLE)",length(data@doseGrid)),
rep("Expected Efficacy",length(data@doseGrid)),
rep("Gain",length(data@doseGrid))),
Type=factor("Estimate",levels="Estimate")
))
plot1 <- ggplot(data=gdata, aes(x=x,y=y))+geom_line(aes(group=group,color=group),size=1.5)+
ggplot2::scale_colour_manual(name="curves",values=c("blue","green3","red"))+
xlab("Dose Level")+ xlim(c(0,max(data@doseGrid)))+
ylab(paste("Values")) + ylim(c(min(gdata$y),max(gdata$y)))
if ((signif(TDtargetEndOfTrialEstimate,4) < LowestDose) |(signif(TDtargetEndOfTrialEstimate,4) > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",EndOfTrialtargetprob*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <-plot1 + geom_point(data=data.frame(x=TDtargetEndOfTrialEstimate,y=EndOfTrialtargetprob),aes(x=x,y=y),colour="violet", shape=16, size=8) +
annotate("text",label=paste(paste("TD",EndOfTrialtargetprob*100),"Estimate"),x=TDtargetEndOfTrialEstimate-3,y=0.2,size=5,colour="violet")}
if ((signif(Gstar,4) < LowestDose) |(signif(Gstar,4) > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste("Estimated Gstar=",paste(Gstar," not within dose Grid")))} else {plot1 <- plot1}
} else {plot1 <- plot1 +
geom_point(data=data.frame(x=Gstar,y=MaxGain),aes(x=x,y=y),colour="green3", shape=17, size=8) +
annotate("text",label="Max Gain Estimate",x=Gstar,y=MaxGain-0.1,size=5,colour="green3")}
mylabel=format(DuringTrialtargetprob,digits=2)
if ((signif(TDtargetDuringTrialEstimate,4) < LowestDose) |(signif(TDtargetDuringTrialEstimate,4) > max(data@doseGrid))) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",DuringTrialtargetprob*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_point(data=data.frame(x=signif(TDtargetDuringTrialEstimate,4),y=DuringTrialtargetprob),aes(x=x,y=y),colour="orange", shape=15, size=8) +
annotate("text",label=paste(paste("TD",DuringTrialtargetprob*100),"Estimate"),x=TDtargetDuringTrialEstimate+25,
y=DuringTrialtargetprob+0.01,size=5,colour="orange")
}
if (doselimit > max(data@doseGrid)) {maxdoselimit <- max(data@doseGrid)} else {maxdoselimit<-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="brown", lwd=1.1) +
annotate("text",label="Max",x=maxdoselimit-2,y=max(gdata$y),size=5,colour="brown")
plot1 <-plot1 +
geom_vline(xintercept=ret, colour="purple", lwd=1.1) +
annotate("text",label="Next", x=ret+1, y=max(gdata$y)-0.05,size=5,color="purple")
## return next best dose and plot
return(list(nextdose=ret,
DLEDuringTrialtarget=DuringTrialtargetprob,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
TDtargetDuringTrialAtDoseGrid=retD,
DLEEndOfTrialtarget=EndOfTrialtargetprob,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=retE,
GstarEstimate=Gstar,
GstarAtDoseGrid=Gstarret,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
CIGstar=CIGstar,
ratioGstar=ratioGstar,
plot=plot1))
})
## =====================================================================================
## the nextBest method based on DLE and efficacy responses with DLE and efficacy samples
## -------------------------------------------------------------------------- ----------
##' @describeIn nextBest for slots \code{nextBest},\code{doselimit}, \code{data} and \code{SIM}. This is
##' a function to find the next best dose based on the 'NextBestMaxGainSamples'
##' class rule. This a method based on the DLE responses and efficacy responses with DLE and
##' efficacy samples. Effmodel must be of class \code{\linkS4class{Effloglog}} or
##' \code{\linkS4class{EffFlexi}}.
##' @param Effsamples the efficacy samples of \code{\linkS4class{Samples}} class object
##'
##' @importFrom ggplot2 ggplot geom_histogram xlab ylab xlim aes geom_vline
##' geom_text
##' @example examples/Rules-method-nextbest_MaxGainSamples.R
setMethod("nextBest",
signature=
signature(nextBest="NextBestMaxGainSamples",
doselimit="numeric",
samples="Samples",
model="ModelTox",
data="DataDual"),
def=
function(nextBest,doselimit,samples,model,data,Effmodel,Effsamples,SIM=FALSE,...){
if(is(Effmodel, "Effloglog"))
{
##first get the probDLE samples
points <- data@doseGrid
probDLESamples <- matrix(nrow=sampleSize(samples@options),
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:
probDLESamples[, i] <- prob(dose=points[i],
model=model,
samples=samples)
}
probDLE <- apply(probDLESamples,2,FUN=nextBest@TDderive)
DuringTrialtargetprob <- nextBest@DLEDuringTrialtarget
EndOfTrialtargetprob <- nextBest@DLEEndOfTrialtarget
## Find the TDtarget samples for During Trial and End of trial
TDtargetEndOfTrialSamples <- dose(prob=EndOfTrialtargetprob,
model=model,
samples=samples)
TDtargetDuringTrialSamples<-dose(prob=DuringTrialtargetprob,
model=model,
samples=samples)
## Find the TDtarget Estimate for During Trial and End of trial
TDtargetEndOfTrialEstimate <- as.numeric(nextBest@TDderive(TDtargetEndOfTrialSamples))
## Ensure the estimate is within dose range
#TDtargetEndOfTrialEstimate <- min(TDtargetEndOfTrialEstimate,max(data@doseGrid))
TDtargetDuringTrialEstimate<-as.numeric(nextBest@TDderive(TDtargetDuringTrialSamples))
## Ensure the estimate is within dose range
#TDtargetDuringTrialEstimate <- min(TDtargetDuringTrialEstimate,max(data@doseGrid))
##we have to get samples from the gain values at all dose levels
ExpEffSamples <- matrix(nrow=sampleSize(Effsamples@options),
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:
ExpEffSamples[, i] <- ExpEff(dose=points[i],
model=Effmodel,
samples=Effsamples)
}
ExpEff <- apply(ExpEffSamples,2,FUN=nextBest@Gstarderive)
GainSamples <- matrix(nrow=sampleSize(samples@options),
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],
DLEmodel=model,
DLEsamples=samples,
Effmodel=Effmodel,
Effsamples=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]
##Obtain the Gstar estimate which is the nth percentile of the Gstar samples
Gstar <- as.numeric(nextBest@Gstarderive(GstarSamples))
##Ensure the estimate is within dose range
#Gstar <- min(Gstar,max(data@doseGrid))
gainvalues <- apply(GainSamples,2,FUN=nextBest@Gstarderive)
dosesOK <- which(data@doseGrid <= doselimit)
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
##FIND the next dose which is the minimum between TDtargetDuringTrial and Gstar
nextdose<-min(TDtargetDuringTrialEstimate,Gstar)
##Find the dose level in doseGrid closest below nextdose
index <- suppressWarnings(max(which((signif(nextdose,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level in doseGrid closest below TDtargetEndOfTrial
indexE <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retE <- data@doseGrid[indexE]
##Find the dose level in doseGrid closest below TDtargetDuringTrial
indexD <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retD <- data@doseGrid[indexD]
##Find the dose level in doseGrid closest below Gstar
Gstarindex <- suppressWarnings(max(which((signif(Gstar,digits=4) - data@doseGrid[dosesOK]) >= 0)))
Gstarret <- data@doseGrid[Gstarindex]
##Find the 95% credibility interval of Gstar and its ratio of the upper to the lower limit
CIGstar <- as.numeric(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 <- as.numeric(quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975)))
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData<-data.frame(dose=rep(data@doseGrid,3),
values=c(probDLE,
ExpEff,
gainvalues))
gdata<-with(plotData,
data.frame(x=dose,
y=values,
group=c(rep("p(DLE)",length(data@doseGrid)),
rep("Expected Efficacy",length(data@doseGrid)),
rep("Gain",length(data@doseGrid))),
Type=factor("Estimate",levels="Estimate")
))
## produce plot
plot1 <- ggplot() +
geom_histogram(data=
data.frame(x=GstarSamples),
aes(x=x),
fill = "darkgreen", colour = "green3", binwidth = 25) +
xlab("Gstar")+ xlim(c(0,max(data@doseGrid)))+
ylab("Posterior density")
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
if (signif(TDtargetDuringTrialEstimate,4) < LowestDose |signif(TDtargetDuringTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",DuringTrialtargetprob*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetDuringTrialEstimate, colour="orange", lwd=1.1) +
annotate("text",label=paste(paste("TD",DuringTrialtargetprob*100),"Estimate"),
x=TDtargetDuringTrialEstimate,y=0,hjust=-0.1, vjust = -20,size=5,colour="orange")}
if (signif(TDtargetEndOfTrialEstimate,4) < LowestDose |signif(TDtargetEndOfTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",EndOfTrialtargetprob*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetEndOfTrialEstimate, colour="violet", lwd=1.1) +
annotate("text",label=paste(paste("TD",EndOfTrialtargetprob*100),"Estimate"),
x=TDtargetEndOfTrialEstimate,y=0,hjust=-0.1, vjust = -25,size=5,colour="violet")}
if (signif(Gstar,4) < LowestDose |signif(Gstar,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste("Estimated Gstar=",paste(Gstar," not within dose Grid")))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=Gstar, colour="green", lwd=1.1) +
annotate("text",label=" Gstar Estimate",
x=Gstar,y=0,hjust=-0.1, vjust = -25,size=5,colour="green")}
if (doselimit > max(data@doseGrid)){maxdoselimit<-max(data@doseGrid)} else {maxdoselimit <-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -35),
colour="red")
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0.1, vjust = -30),
colour="blue")
## return next best dose and plot
return(list(nextdose=ret,
DLEDuringTrialtarget=DuringTrialtargetprob,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
TDtargetDuringTrialAtDoseGrid=retD,
DLEEndOfTrialtarget=EndOfTrialtargetprob,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=retE,
GstarEstimate=Gstar,
GstarAtDoseGrid=Gstarret,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
CIGstar=CIGstar,
ratioGstar=ratioGstar,
plot=plot1))
} else if(is(Effmodel, "EffFlexi")) {
##first get the probDLE samples
points <- data@doseGrid
probDLESamples <- matrix(nrow=sampleSize(samples@options),
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:
probDLESamples[, i] <- prob(dose=points[i],
model=model,
samples=samples)
}
probDLE <- apply(probDLESamples,2,FUN=nextBest@TDderive)
DuringTrialtargetprob <- nextBest@DLEDuringTrialtarget
EndOfTrialtargetprob <- nextBest@DLEEndOfTrialtarget
## Find the TDtarget samples for During Trial and End of trial
TDtargetEndOfTrialSamples <- dose(prob=EndOfTrialtargetprob,
model=model,
samples=samples)
TDtargetDuringTrialSamples<-dose(prob=DuringTrialtargetprob,
model=model,
samples=samples)
## Find the TDtarget Estimate for During Trial and End of trial
TDtargetEndOfTrialEstimate <- as.numeric(nextBest@TDderive(TDtargetEndOfTrialSamples))
## Ensure the estimate is within dose range
#TDtargetEndOfTrialEstimate <- min(TDtargetEndOfTrialEstimate,max(data@doseGrid))
TDtargetDuringTrialEstimate<-as.numeric(nextBest@TDderive(TDtargetDuringTrialSamples))
## Ensure the estimate is within dose range
#TDtargetDuringTrialEstimate <- min(TDtargetDuringTrialEstimate,max(data@doseGrid))
##we have to get samples from the gain values at all dose levels
ExpEffsamples <- Effsamples@data$ExpEff
ExpEff <- apply(ExpEffsamples,2,FUN=nextBest@Gstarderive)
GainSamples <- matrix(nrow=sampleSize(samples@options),
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],
DLEmodel=model,
DLEsamples=samples,
Effmodel=Effmodel,
Effsamples=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]
##Obtain the Gstar estimate which is the 50th percentile of the Gstar samples
Gstar <- as.numeric(nextBest@Gstarderive(GstarSamples))
##Ensure the estimate is within dose range
#Gstar <- min(Gstar,max(data@doseGrid))
gainvalues <- apply(GainSamples,2,FUN=nextBest@Gstarderive)
dosesOK <- which(data@doseGrid <= doselimit)
## For placebo design, if safety allow, exclude placebo from
## the recommended next doses
if(data@placebo & (length(dosesOK) > 1L) ){
dosesOK <- dosesOK[-1]
}
##FIND the next dose which is the minimum between TDtargetDuringTrial and Gstar
nextdose<-min(TDtargetDuringTrialEstimate,Gstar)
##Find the dose level in doseGrid closest below nextdose
index <- suppressWarnings(max(which((signif(nextdose,digits=4) - data@doseGrid[dosesOK]) >= 0)))
ret <- data@doseGrid[dosesOK][index]
##Find the dose level in doseGrid closest below TDtargetEndOfTrial
indexE <- suppressWarnings(max(which((signif(TDtargetEndOfTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retE <- data@doseGrid[indexE]
##Find the dose level in doseGrid closest below TDtargetDuringTrial
indexD <- suppressWarnings(max(which((signif(TDtargetDuringTrialEstimate,digits=4) - data@doseGrid[dosesOK]) >= 0)))
retD <- data@doseGrid[indexD]
##Find the dose level in doseGrid closest below Gstar
Gstarindex <- suppressWarnings(max(which((signif(Gstar,digits=4) - data@doseGrid[dosesOK]) >= 0)))
Gstarret <- data@doseGrid[Gstarindex]
##Find the 95% credibility interval of Gstar and its ratio of the upper to the lower limit
CIGstar <- as.numeric(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 <- as.numeric(quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975)))
ratioTDEOT <- as.numeric(CITDEOT[2]/CITDEOT[1])
plotData<-data.frame(dose=rep(data@doseGrid,3),
values=c(probDLE,
ExpEff,
gainvalues))
gdata<-with(plotData,
data.frame(x=dose,
y=values,
group=c(rep("p(DLE)",length(data@doseGrid)),
rep("Expected Efficacy",length(data@doseGrid)),
rep("Gain",length(data@doseGrid))),
Type=factor("Estimate",levels="Estimate")
))
## produce plot
plot1 <- ggplot() +
geom_histogram(data=
data.frame(x=GstarSamples),
aes(x=x),
fill = "darkgreen", colour = "green3", binwidth = 25) +
xlab("Gstar")+ xlim(c(0,max(data@doseGrid)))+
ylab("Posterior density")
#if(data@placebo) {
#n <- length(data@doseGrid)
#LowestDose <- sort(data@doseGrid)[2]} else {
LowestDose <- min(data@doseGrid)
#}
if (signif(TDtargetDuringTrialEstimate,4) < LowestDose |signif(TDtargetDuringTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",DuringTrialtargetprob*100),paste("=",paste(TDtargetDuringTrialEstimate," not within dose Grid"))))
} else {plo1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetDuringTrialEstimate, colour="orange", lwd=1.1) +
annotate("text",label=paste(paste("TD",DuringTrialtargetprob*100),"Estimate"),
x=TDtargetDuringTrialEstimate,y=0,hjust=-0.1, vjust = -20,size=5,colour="orange")}
if (signif(TDtargetEndOfTrialEstimate,4) < LowestDose |signif(TDtargetEndOfTrialEstimate,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste(paste("Estimated TD",EndOfTrialtargetprob*100),paste("=",paste(TDtargetEndOfTrialEstimate," not within dose Grid"))))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=TDtargetEndOfTrialEstimate, colour="violet", lwd=1.1) +
annotate("text",label=paste(paste("TD",EndOfTrialtargetprob*100),"Estimate"),
x=TDtargetEndOfTrialEstimate,y=0,hjust=-0.1, vjust = -25,size=5,colour="violet")}
if (signif(Gstar,4) < LowestDose |signif(Gstar,4) > max(data@doseGrid)) {
if (SIM==FALSE){
plot1<-plot1
print(paste("Estimated Gstar=",paste(Gstar," not within dose Grid")))
} else {plot1 <- plot1}
} else {
plot1 <- plot1+
geom_vline(xintercept=Gstar, colour="green", lwd=1.1) +
annotate("text",label=" Gstar Estimate",
x=Gstar,y=0,hjust=+0.6, vjust = -25,size=5,colour="green")}
if (doselimit > max(data@doseGrid)){maxdoselimit<-max(data@doseGrid)} else {maxdoselimit <-doselimit}
plot1 <- plot1 +
geom_vline(xintercept=maxdoselimit, colour="red", lwd=1.1) +
geom_text(data=
data.frame(x=maxdoselimit),
aes(x, 0,
label = "Max", hjust = +1, vjust = -35),
colour="red")
plot1 <- plot1 +
geom_vline(xintercept=ret, colour="blue", lwd=1.1) +
geom_text(data=
data.frame(x=ret),
aes(x, 0,
label = "Next", hjust = 0.1, vjust = -30),
colour="blue")
## return next best dose and plot
return(list(nextdose=ret,
DLEDuringTrialtarget=DuringTrialtargetprob,
TDtargetDuringTrialEstimate=TDtargetDuringTrialEstimate,
TDtargetDuringTrialAtDoseGrid=retD,
DLEEndOfTrialtarget=EndOfTrialtargetprob,
TDtargetEndOfTrialEstimate=TDtargetEndOfTrialEstimate,
TDtargetEndOfTrialAtDoseGrid=retE,
GstarEstimate=Gstar,
GstarAtDoseGrid=Gstarret,
CITDEOT=CITDEOT,
ratioTDEOT=ratioTDEOT,
CIGstar=CIGstar,
ratioGstar=ratioGstar,
plot=plot1))
} else stop("Effmodel needs to be of class Effloglog or EffFlexi")
})
## ------------------------------------------------------------------------------------------------
## 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("stopTrial",
signature=
signature(stopping="StoppingTDCIRatio",
dose="ANY",
samples="Samples",
model="ModelTox",
data="ANY"),
def=
function(stopping,dose,samples,model,data,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##check id targetEndOfTrial is a probability
stopifnot(is.probability(targetEndOfTrial))
## find the TDtarget End of Trial samples
TDtargetEndOfTrialSamples <- dose(prob=targetEndOfTrial,
model=model,
samples=samples)
##Find the upper and lower limit of the 95% credibility interval
CI <- quantile(TDtargetEndOfTrialSamples, probs=c(0.025,0.975))
##The ratio of the upper to the lower 95% credibility interval
ratio <- as.numeric(CI[2]/CI[1])
##so can we stop?
doStop <- ratio <= stopping@targetRatio
##generate messgae
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"),
"targetRatio =", stopping@targetRatio)
##return both
return(structure(doStop,
messgae=text))
})
## ----------------------------------------------------------------------------------------------
## 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,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##check if targetEndOfTrial is a probability
stopifnot(is.probability(targetEndOfTrial))
## find the TDtarget End of Trial
TDtargetEndOfTrial <- dose(prob=targetEndOfTrial,
model=model)
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(targetEndOfTrial/(1-targetEndOfTrial))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval
CI <- c()
CI[2] <- exp(log(TDtargetEndOfTrial) + 1.96* sqrt(varEta))
CI[1] <- exp(log(TDtargetEndOfTrial) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratio <- as.numeric(CI[2]/CI[1])
##so can we stop?
doStop <- ratio <= stopping@targetRatio
##generate messgae
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"),
"targetRatio =", stopping@targetRatio)
##return both
return(structure(doStop,
messgae=text))
})
## --------------------------------------------------------------------------------------------------
## 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="StoppingGstarCIRatio",
dose="ANY",
samples="Samples",
model="ModelTox",
data="DataDual"),
def=
function(stopping,dose,samples,model,data,TDderive, Effmodel,Effsamples,Gstarderive,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##checks
stopifnot(is.probability(targetEndOfTrial))
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(prob=targetEndOfTrial,
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=sampleSize(samples@options),
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],
DLEmodel=model,
DLEsamples=samples,
Effmodel=Effmodel,
Effsamples=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@targetRatio
##generate messgae
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,"TDatrgetEndOfTrial estimate","Gstar estimate"), "is smaller with ratio =",
round(ratio,4), " which is " , ifelse(doStop,"is less than or equal to","greater than"),
"targetRatio =", stopping@targetRatio)
text <- c(text1,text2,text3)
##return both
return(structure(doStop,
message=text))
})
## -----------------------------------------------------------------------------------------------
## 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="StoppingGstarCIRatio",
dose="ANY",
samples="missing",
model="ModelTox",
data="DataDual"),
def=
function(stopping,dose,model,data,Effmodel,...){
targetEndOfTrial <- stopping@targetEndOfTrial
##checks
stopifnot(is.probability(targetEndOfTrial))
stopifnot(is(Effmodel,"ModelEff"))
## find the TDtarget End of Trial
TDtargetEndOfTrial <- dose(prob=targetEndOfTrial,
model=model)
##Find the dose with maximum gain value
Gainfun<-function(DOSE){
-gain(DOSE,DLEmodel=model,Effmodel=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@c)} 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 <- t(deltaG)%*%covBETA%*%deltaG
##Find the upper and lower limit of the 95% credibility interval of Gstar
CIGstar <-c()
CIGstar[2] <- exp(logGstar + 1.96* sqrt(varlogGstar))
CIGstar[1] <- exp(logGstar - 1.96* sqrt(varlogGstar))
##The ratio of the upper to the lower 95% credibility interval
ratioGstar <- as.numeric(CIGstar[2]/CIGstar[1])
##Find the variance of the log of the TDtargetEndOfTrial(eta)
M1 <- matrix(c(-1/(model@phi2), - (log(targetEndOfTrial/(1-targetEndOfTrial))-model@phi1)/(model@phi2)^2),1,2)
M2 <- model@Pcov
varEta <- M1%*%M2%*%t(M1)
##Find the upper and lower limit of the 95% credibility interval of
##TDtargetEndOfTrial
CITDEOT <- c()
CITDEOT[2] <- exp(log(TDtargetEndOfTrial) + 1.96* sqrt(varEta))
CITDEOT[1] <- exp(log(TDtargetEndOfTrial) - 1.96* sqrt(varEta))
##The ratio of the upper to the lower 95% credibility interval
ratioTDEOT <- as.numeric(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@targetRatio
##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"),
"targetRatio =", stopping@targetRatio)
text <- c(text1,text2,text3)
##return both
return(structure(doStop,
message=text))
})
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.