In 0liver0815/onc-crmpack-test: Object-Oriented Implementation of CRM Designs
Date r params$d
knitr::opts_chunk$set(echo = TRUE,cache =FALSE)
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols: Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
###
library(vcd)
library(MASS)
library("eha")
library("msm")
library("crmPack")
source("safetywindow.r")
source("mDA-CRM-model.r")
source("mDA-CRM-design.r")
library (xlsx)
##;
############ rules agreed with team
## Target toxicity interval: 20% – 35% DLT rate.
## Starting dose and dose range: 1.6 mg (to be confirmed) up to 400 mg (3 weekly flat dosing)
## Cohort size : at least 3 patients per cohort
## Maximum dose increments
## 100% increase prior to 1st DLT
## 50% increase after 1st DLT
## <25% probability to overdose (to be above targeted toxicity interval)
## Stopping criteria
## Chance (or likelihood) of being in the targeted toxicity interval >40%
## 6 subjects or at least 2 cohorts having dose within 20% dose of the MTD
## Clinical judgment: no benefit with additional cohorts
## A maximum of 60 patients is reached
## DLT window and safety windows
## DLT window 28 days from first administration.
## 1st patient safety window 1 week before enrolling the remaining patients,
## 1 day interval between 2nd and 3rd patient.
## for 2nd and 3rd patient at least 2 weeks safety window before starting the next cohort.
## Intra-patient escalation: allowed only if part of a dosing regimen explored, otherwise not planned
## from ceatcb safe dose it has been applied a 20 fold safety factor for determining the starting dose for cam5
## starting dose for single patient cohort part is 65 microg
## doses of cea tcb are five fold higher than estimated dose of ceacam5 tcb so need to adjust to
## reflect the difference in potency
myTmax <-28 #Tmax is the lenght of the DLT period in this project (CAM5)
mynpiece <- c(myTmax/7)
myStartDose <- 1.6
###
mydoseGrid=c(seq(from=0.5, to=10, by=0.1),
seq(from=11, to=20, by=1),25,27,
seq(from=30, to=100, by=5),120,
seq(from=150, to=400, by=50))
set.seed(23)
################## Prior data from CEATCB ########################
## Monotherapy study, N = 80 patients received at least one dose of IV RO6958688 QW up to 600 mg in
## escalation part of the trial (MTD for C1/C2 QW defined as 400 mg). Safety profile dominated by
## IRRs and Tumor Flare.
## (a dose grid up to 80 mg with 1 mg spacing, and from 85 to 1000 mg with 5 mg spacing,will be used).
## Five DLTs in five patients, at doses over 40 mg occurring shortly after the first administration
## 2 DLTs @ 600 mg (N = 2 pts), Respiratory failure Gr5, Colitis Gr4 (discontinued)
## 1 DLT @ 300 mg (N = 8 pts), Diarrhea Gr3 (resolved, C1D5-C1D18)
## 1 DLT @ 60 mg (N = 7 pts), Hypoxia Gr3 (resolved, C1D2-C1D10)
## 1 DLT @ 40mg (N = 12 pts), Dyspnea Gr3 (resolved, C1D2-C1D7)
## from protocol prior parameters for tcb model without pretreatment
## doseGrid 0.05 1 50 100 500 1000
## mean 1.021, -0.682
## cov 1.887,0.023,0.023,0.001
## reference dose 500
#### reading in the tcb dataset
dataTCB<- read.xlsx(file="dltdataTCB.xlsx",1)
dataTCB$DLTDY[dataTCB$DLTDY>myTmax] <- myTmax ## note DLTDY cannot be bigger than the dlt window used later in the model
dataTCB$DLTDY[dataTCB$DLTDY==0] <- 0.1
#looking at original set up for cea tcb
mydoseGrid0 <- c(2.5, seq (3,85,1),seq(90,1000,5)) #from protocol
myrefDose0<- c(500)
model0 <- crmPack:::LogisticLogNormal(mean=c(1.021, -0.682),
cov= matrix(c(1.887 , 0.023, 0.023, 0.001), nrow=2),
refDose=myrefDose0) # keeping on the scale that is needed in this project
#reading the data from ceatcb ;
data0 <- crmPack:::Data(x=c(dataTCB$EXDOSE),y=c(dataTCB$DLT), doseGrid=mydoseGrid0)
options <- crmPack:::McmcOptions(burnin=10000,
step=2,
samples=50000)
set.seed(94)
samples0 <- crmPack:::mcmc (data0,model0,options)
plot(samples0,model0,data0)
### WARNING THE MTD IN HEREIS NOT AROUND 400.mtc was edfine also on the basis of other clinical consideration and the model prior was too informative
####
# Approach 1;
# Based on the mixture prior to approximate a lognormal prior;
### getting mean and cov priors from mixture of informative and non informative components for ceatcb:
myrefDose<- c(500/5) #scaling ref dose as well the doses from tcb dataset
model_TCBS <- crmPack:::LogisticLogNormal(mean=c(1.021, -0.682),
cov= matrix(c(1.887 , 0.023, 0.023, 0.001), nrow=2),
refDose=myrefDose) # keeping on the scale that is needed in this project
#reading the data from ceatcb scaling the doses;
data_TCBS <- crmPack:::Data(x=c(dataTCB$EXDOSE/5),y=c(dataTCB$DLT), doseGrid=mydoseGrid)
options <- crmPack:::McmcOptions(burnin=10000,
step=2,
samples=50000)
set.seed(94)
samples_TCBS <- crmPack:::mcmc (data_TCBS,model_TCBS,options)
plot(samples_TCBS,model_TCBS,data_TCBS)
```{R inf_mod}
informative component
inmodel<- approximate(samples_TCBS, model_TCBS, data_TCBS, points = c(1,20,50,100,200,400),
refDose = myrefDose, logNormal = TRUE, verbose = TRUE)
creating an empty dataset, provide only the dose grid:
emptydata_inm <- Data(doseGrid=mydoseGrid)
obtain prior samples with this Data object note inmodel is with posterior parameters !
priorsamples_inm <- mcmc(emptydata_inm, inmodel, options)
## then produce the plot
plot(priorsamples_inm, inmodel, emptydata_inm)
pi<- plot(priorsamples_inm, inmodel, emptydata_inm) #saving the plot as object to plot it later with multiplot
```r
## minimal inf component
coarseGrid <- c(1,20,50,100,200,400)
minInfModel <- MinimalInformative(dosegrid = coarseGrid,
refDose=myrefDose,
threshmin=0.1,
threshmax=0.3,
logNormal = TRUE)
minInfModel$model@cov
## be carefyll that these two plots are one on the top of the other so check specification for output!
matplot(x=coarseGrid,
y=minInfModel$required,
type="b", pch=19, col="blue", lty=1,
xlab="dose",
ylab="prior probability of DLT")
matlines(x=coarseGrid,
y=minInfModel$quantiles,
type="b", pch=19, col="red", lty=1)
legend("right",
legend=c("quantiles", "approximation"),
col=c("blue", "red"),
lty=1,
bty="n")
noninmodel<-minInfModel$model
emptydata <- Data(doseGrid=mydoseGrid)
priorsamples_min <- mcmc(emptydata, noninmodel, options)
## then produce the plot
plot(priorsamples_min, noninmodel, emptydata)
pni<-plot(priorsamples_min, noninmodel, emptydata)#saving the plot as object to plot it later with multiplot
## mix informative and minimal informative models exploring different weights;
components<- list(inf=list(mean=inmodel@mean,cov=inmodel@cov),
noninf=list(mean=noninmodel@mean,cov=noninmodel@cov))
#for weight use the same order for the model as above, exploring different weights;
mixmodel_1090 <- LogisticNormalFixedMixture(components=components, weights=c(10,90), refDose=myrefDose, logNormal = TRUE)
mixmodel_2575 <- LogisticNormalFixedMixture(components=components, weights=c(25,75), refDose=myrefDose, logNormal = TRUE)
mixmodel_4060 <- LogisticNormalFixedMixture(components=components, weights=c(40,60), refDose=myrefDose, logNormal = TRUE)
mixmodel_7030 <- LogisticNormalFixedMixture(components=components, weights=c(70,30), refDose=myrefDose, logNormal = TRUE)
emptydata <- Data(doseGrid=mydoseGrid)
priorsamples_1090 <- mcmc(emptydata, mixmodel_1090, options)
plot(priorsamples_1090, mixmodel_1090, emptydata)
## then produce the plot
p1<- plot(priorsamples_1090, mixmodel_1090, emptydata)
priorsamples_2575 <- mcmc(emptydata, mixmodel_2575, options)
plot(priorsamples_2575, mixmodel_2575, emptydata)
## then produce the plot
p2<- plot(priorsamples_2575, mixmodel_2575, emptydata)
priorsamples_4060 <- mcmc(emptydata, mixmodel_4060, options)
plot(priorsamples_4060, mixmodel_4060, emptydata)
## then produce the plot
p3<- plot(priorsamples_4060, mixmodel_4060, emptydata)
priorsamples_7030 <- mcmc(emptydata, mixmodel_7030, options)
plot(priorsamples_7030, mixmodel_7030, emptydata)
## then produce the plot
p4<- plot(priorsamples_7030, mixmodel_7030, emptydata)
#more plots on one page
multiplot(pi, pni, p1, cols=1)
multiplot(p1, p2, p3, cols=1)
multiplot(p1, p3, p4, cols=1)
#overlaying plots
p1innoin <- p1 +
geom_line(data=pni$data,
aes(x=x, y=y, group=group, linetype=Type),
colour="blue"
)+
geom_line(data=pi$data,
aes(x=x, y=y, group=group, linetype=Type),
colour="black")
#make sure the caption is correct!!
p1innoin+ labs (caption =c("In black the Inf. model, in blue the min. inf. mod. and in red the mixture with weight 10/90") )
p1234 <- p1 +
geom_line(data=p2$data,
aes(x=x, y=y, group=group, linetype=Type),
colour="blue")+
geom_line(data=p3$data,
aes(x=x, y=y, group=group, linetype=Type),
colour="purple")+
geom_line(data=p4$data,
aes(x=x, y=y, group=group, linetype=Type),
colour="black")
p1234 +labs (caption =c("Mixed models weights (inf/min.inf.): in red 10/90, in blue 25/75, in purple 40/60, in black 70/30") )
## need to decide on a mixture prior for next steps and get the parameters
## get the 25/75
#overlaying plots
p2innoin <- p2 +
geom_line(data=pni$data,
aes(x=x, y=y, group=group, linetype=Type),
colour="blue"
)+
geom_line(data=pi$data,
aes(x=x, y=y, group=group, linetype=Type),
colour="black")
#make sure the caption is correct!!
p2innoin+ labs (caption =c("In black the Inf. model, in blue the min. inf. mod. and in red the mixture with weight 25/75") )
posterior_hist <- approximate(priorsamples_2575, mixmodel_2575, emptydata, refDose=mixmodel_2575@refDose, logNormal = TRUE,
seed = 12345,
control = list(threshold.stop = 0.01, maxit = 50000,
temperature = 50000, max.time = 120))
posterior_hist@mean
posterior_hist@cov
posterior_hist@refDose
emptydata <- Data(doseGrid=mydoseGrid)
priorsamples_post <- mcmc(emptydata, posterior_hist, options)
plot_priorsam <- plot(priorsamples_post , posterior_hist, emptydata)
multiplot(p2,plot_priorsam)
#posterior_hist2 <- approximate(priorsamples_2575, mixmodel_2575, emptydata, #refDose=mixmodel_2575@refDose, logNormal = FALSE)
#posterior_hist2@mean
#posterior_hist2@cov
#posterior_hist2@refDose
#emptydata <- Data(doseGrid=mydoseGrid)
#priorsamples_post2 <- mcmc(emptydata, posterior_hist2, options)
#plot_priorsam2 <- plot(priorsamples_post2 , posterior_hist2, emptydata)
# multiplot(p2,plot_priorsam,plot_priorsam2)
####Approach 2;
####Use quantiles to log normal function to estimate the mixture prior;
# test <- Quantiles2LogisticNormal(dosegrid=c(1, 10, 30, 40, 80),
# refDose=56,
# lower= c(0.01,0.02,0.03,0.05,0.12),
# median= c(0.02,0.08,0.3,0.45,0.77), #median probability of DLT at each dose in dosegrid is enough to put 5 to 6 doses,
## these are derived from the posterior logistic dose
# upper= c(0.10,0.6,0.92,0.97,0.99),
# logNormal=TRUE)
######################### modelling both dose and time ###############################
## reading in data including time of DLT MAKE SURE there are no 0 time!!
data <- DataDA(x=c(dataTCB$EXDOSE)/5, ### patient dose level doses of cea tcb are five fold higher than estimated dose of ceacam5 tcb so need to adjust
y=c(dataTCB$DLT), ### patient DLT 1=yes 0 =no
doseGrid=mydoseGrid,
u=c(dataTCB$DLTDY), ###u= DLT free survival values time to dlt or time to censoring if no dlt.
## NOTE u entries must not be larger than Tmax and
### dlt occuring outside the considered dlt window should be excluded.
t0=c(dataTCB$dayd), ### time from recruitment of first patient to recruitment of each of the remaining patients
Tmax=myTmax,
npiece=mynpiece)
#2) Structure of the model class
#mean=c(-0.85,1)
#cov=matrix(c(1,-0.5,-0.5,1),nrow=2)
#refDose=56
#non-informative prior for lambdas (number of pieces=npiece_)
#l=as.numeric(t(apply(as.matrix(c(1:npiece_),1,npiece_),2,lambda_prior)))
##need to fill in
npiece_=mynpiece
Tmax_=myTmax
## note dose grid in empty data need to be same as above
emptydata_mda <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax,npiece=mynpiece)
### do not change the lambda prior is based on paper (is a non non informative lambda)
lambda_prior<-function(k){
npiece_/(Tmax_*(npiece_-k+0.5))
}
### for lambda if want to use informative then may be better to fit from exponential otherwise may be the posterior of lambda is not stable enough
### coud try both with informative and with noninf lambda and see the model performance in the simulations
##try non informative first
mymda_model<-DALogisticLogNormal(mean=c(posterior_hist@mean[1],posterior_hist@mean[2]), ##priors parameters use posterior derived from tcb mixture
cov=matrix(c(posterior_hist@cov[1],posterior_hist@cov[2],posterior_hist@cov[3],posterior_hist@cov[4]),nrow=2),
refDose=posterior_hist@refDose,
l=as.numeric(t(apply(as.matrix(c(1:npiece_),1,npiece_),2,lambda_prior))),
##lambda here not informative and assuming higher prob of event later in the dlt time
C_par=2) ## do not change C_par =2 it control the credible intervals around the H rate.
### ##priors parameters use posterior derived from tcb mixture
#TRY DIFFRENT LAMBDA TO SEE THE EFFECT
mymda_model2<-DALogisticLogNormal(mean=c(posterior_hist@mean[1],posterior_hist@mean[2]), cov=matrix(c(posterior_hist@cov[1],posterior_hist@cov[2],posterior_hist@cov[3],posterior_hist@cov[4]),nrow=2),
refDose=posterior_hist@refDose,
#l=as.numeric(t(apply(as.matrix(c(1:npiece_),1,npiece_),2,lambda_prior))),
l=rev(c(0.04081633, 0.05714286, 0.09523810, 0.28571429)),
##lambda here assuming higher probability at the beginning of the dlt period
C_par=2) ## do not change C_par =2 it control the credible intervals around the H rate.
options <- crmPack:::McmcOptions(burnin=10000,
step=2,
samples=50000)
#3) Obtain the posterior using existing data for example
#set.seed(94)
#samples <- mcmc (data,mymda_model,options)
#4) use ggmcmc to diagnose
library(ggmcmc)
#alpha0samples=get(samples,"alpha0")
#print(ggs_traceplot(alpha0samples))
#print(ggs_autocorrelation(alpha0samples))
#5) plot the model fit
#plot(samples, mymda_model,data,hazard=TRUE)
# plot(samples, model,data,hazard=FALSE)#option FALSE does not work!
## PLOTTING PRIOR mDA
emptydata <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax, npiece=mynpiece)
Priorsamples <- mcmc (emptydata,mymda_model,options)
alpha0samples=get(Priorsamples,"alpha0")
print(ggs_traceplot(alpha0samples))
print(ggs_autocorrelation(alpha0samples))
#plot the model fit
plot(Priorsamples, mymda_model,emptydata,hazard=TRUE)
plotPrmda <- plot(Priorsamples, mymda_model,emptydata,hazard=TRUE)
emptydata <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax, npiece=mynpiece)
Priorsamples2 <- mcmc (emptydata,mymda_model2,options)
alpha0samples=get(Priorsamples2,"alpha0")
print(ggs_traceplot(alpha0samples))
print(ggs_autocorrelation(alpha0samples))
#plot the model fit
plot(Priorsamples2, mymda_model2,emptydata,hazard=TRUE)
plotPrmda2 <- plot(Priorsamples2, mymda_model2,emptydata,hazard=TRUE)
#6) Escalation rules
##need to fill in (use the same rule in the section 8 of "using the package crmPack: introductory examples")
myIncrements <- IncrementsRelativeDLT(DLTintervals=c(0, 1),
increments=c(1, 0.5))
myNextBest <- NextBestNCRM(target=c(0.2,0.35),
overdose=c(0.35,1),
maxOverdoseProb=0.25)
mySize<- CohortSizeConst(size=3)
myStopping1 <- StoppingTargetProb(target=c(0.2, 0.35),
prob=0.4)
myStopping2 <- StoppingMinPatients(nPatients=60) #should this be 60 instead??
myStopping3 <- StoppingPatientsNearDose(nPatients=6,percentage=20)
myStopping <- (myStopping1 | myStopping2 | myStopping3)
#7) recommended dose for the next cohort given the data
##### but data here are the tcb data..is this what is intended??is just an example of how it will behave?
#nextMaxDose <- maxDose(myIncrements,data=data) # this depend on data!!
#doseRecommendation <- nextBest(myNextBest,
# doselimit=nextMaxDose,
# samples=samples, ##need to run the mcm to get the posterior with current data
# model=mymda_model,
# data=data)
#doseRecommendation$plot
#doseRecommendation$value
##Example 2: run a simulation to evaluate design operating characters;
#1) set up safety window and DADesign
mysafetywindow=SafetyWindowConst(c(7,1),14,14) ####I could use SafetyWindowSize or SafetyWindowConst
### for each cohort I want to follow first patient for 7 days then all for at ### least 14 days
#SafetyWindowConst(patientGap = c(7,3),
# patientFollow = 7, this is the minimal required condition for all patienst in the cohort
# patientFollowMin = 14) min is condition for at least one patient
emptydata <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax, npiece=mynpiece)
design <- DADesign(model=mymda_model,
increments=myIncrements,
nextBest=myNextBest,
stopping=myStopping,
cohortSize=mySize,
data=emptydata,
safetyWindow=mysafetywindow,
startingDose=c(myStartDose))
trial_bh <- examine(design,options=options)
## need to save output of trial_bh and run it overnight is slow!!
## some notes:
## DLTearly_1=1 indicates the scenario that the patients with longest follow up have DLTs, DLTearly_1=2 indicates The patients with shortest follow up have DLTs
#question: it seems that at a certain point trial stop si always = true even in absence of dl and even if the next dose is higher than the current one why?
#answer: In order to review more scenarios, the examine function stops only when newDose <= thisDose or the maximum dose in the dose grid reached. if stop=TRUE that means the stopping rule is fulfilled (this may due to a relatively informative prior). A quick way to check could be: write down the 0 DLT data (till 40mg for example) and run with the mcmc function.
#Also why does it get to up to 6 dlt? are these the sum of the dlt on the current dose and the previous at teh end of the dlt period?
#Yes, sum of the dlt on the current dose and the previous dose who are under the DLT follow up period, i.e. the nfollow I mentioned in the examine code.
# also why even at the initial dose with 3 dlt keep going up?
# This may due to a strong study prior. A way may confirm the reason is to run a nCRM model with the same alpha parameters.
print(proc.time() - t1)
#saving the file (is saving to the working directory)
save(trial_bh,file="trial_bh.RData")
trial_0 <- trial_bh[which (trial_bh$DLTs==0),]
#setting a data set with 0 dlt assuming time of dosing between cohorts is as in safety window
# need a data set with the following variable
# $EXDOSE
# $DLT
# $DLTDY
# $dayd
###################TO DO !!!!#####
#simulation with true scenarios
#2)set up truth curves
myTruth<-function(dose){
model@prob(dose,alpha0=2,alpha1=3)
}
curve(myTruth(x), from=0, to=c(max(mydoseGrid)), ylim=c(0, 1))
myTruth1<-function(dose){
model@prob(dose,alpha0=0.05,alpha1=2)
}
curve(myTruth1(x), from=0, to=c(max(mydoseGrid)), ylim=c(0, 1))
###
onset <- 15 ## median onset of dlt time feed in function below
## need to double check with the value of the exponential sample 1000 and see where is the mean
mytrueTmax<- 28 #Truemax is used to generate the DLT data in the simulation and cannot be lower than onset, so any dlt is contained within 0 and mytrueTmax also feed in the function below. it could be set also bigger than the tmax used in the actual model to see what would happen if the dlt window used is in reality too short.
# if the safety window is set as long as the dlt window then the simulation will evaluate a crm EWOC without mDA so then if the dlt window is shorter than truetmax it is possible to evaluate the advantage of da over other approache so the advantage of using in this case a longer dlt window but with short safety window
exp_cond.cdf<-function(x){
1-(pexp(x,1/onset,lower.tail=FALSE)-pexp(mytrueTmax,1/onset,lower.tail=FALSE))/pexp(mytrueTmax,1/onset)
}
#3) set up simulation settings
mySims <- simulate(design,
args=NULL,
truthTox=myTruth,
truthSurv=exp_cond.cdf,#piece_exp_cond.cdf,
trueTmax=mytrueTmax, ### Truemax is used to generate the DLT data in the simulation. You can image if DLT is generated from
#day 1 to day 28 in the simulation, but tmax=14 days, the design will not capture the full DLT information and it
#leads to underestimation of the MTD.
nsim=10,
seed=819,
mcmcOptions=options,
firstSeparate=TRUE,
deescalate=FALSE, ## it refer to patients already enrolled so with false if a dlt occur in a previous dose or current dose level the patients stay at the level they are.
## if set to TRUE if they deescalate they will be censored at the time of the occurrence of the dlt which cause the decision to de-escalate.
parallel=FALSE) ## NOTE THE OPTION TRUE DOES NOT WORK YET FOR THIS mDA CRM FUNCTION!!!!
# system.time(simulate(design,
# args=NULL,
# truthTox=myTruth,
# truthSurv=exp_cond.cdf,#piece_exp_cond.cdf,
# trueTmax=80,
# nsim=500,
# seed=819,
# mcmcOptions=options,
# firstSeparate=TRUE,
# deescalate=FALSE,
# parallel=FALSE))
#4) interpreate simulation result
#use a similar way as section 9.2 in the "using the package crmPack: introductory examples" document
a=summary(mySims,truth=myTruth)
plot(mySims)
mySims@stopReasons[[2]]
savePlot <- function(myPlot,name) {
png(filename = paste(Sys.Date(),"C:/Users/liaoz4/Documents/R/simulation_result/",name,".png",sep=""), width = 480, height = 480)
print(myPlot)
dev.off()}
###### this part is work on progress!!
myOut<- function (model=mymda_model, Tmax=myTmax, npiece=mynpiece, doseGrid=mydoseGrid){
myout <- data.frame (i=c(1:60),EXDOSE=c(1.6,1.6,1.6,rep(NA,57)),DLT=rep(NA,60),DLTDY=c(rep(NA,60)) ,increment=rep(NA,60), x0DLT=c(1.6,1.6,1.6,rep(NA,57)))
DLTDY=c(22,15,14)
for (i in seq(4,58,3) ) rowSums(dat[,c("b", "c")], na.rm=TRUE)
{
#compute time from 1st patient in
chortind <- (i-1)/3
timeind <- c(rep (DLTDY,each = 1, times = chortind))
myout$DLTDY[1:(i-1)][is.na(myout$DLTDY[1:(i-1)])] <- 0
myout$DLTDY[1:(i-1)] <- myout$DLTDY[1:(i-1)] + timeind[1:(i-1)]
myout$DLTDY[myout$DLTDY>Tmax] <- Tmax #entries must not be larger than Tmax
#computing the dose for the three patient cohort
data <- DataDA(x=c(myout$EXDOSE),
y=c(myout$DLT), ### patient DLT 1=yes 0 =no
doseGrid=doseGrid,
u=c(myout$DLTDY), ###u= DLT free survival values time to dlt or time to censoring if no dlt.
## NOTE u entries must not be larger than Tmax and
### dlt occuring outside the considered dlt window should be excluded.
t0=c(myout$dayd), ### time from recruitment of first patient to recruitment of each of the remaining patients
Tmax=Tmax,
npiece=npiece)
nextBest(myNextBest,
doselimit=nextMaxDose,
samples=samples, ##need to run the mcm to get the posterior with current data
model=model,
data=data)
#myout$EXDOSE[i:(i+2)]<- rep( nextDose (myout$EXDOSE[1:i-1], myout$DLT[1:i-1], model),3)
# to have the table in the column format
myout$x0DLT[(i-3):(i-1)]<-myout$EXDOSE[i:(i+2)]
myout$DLT[1:(i-1)] <- rep (0,c(i-1)) #reset dlt to 0
myout$increment[(i-3):(i-1)] <-rep( round((((myout$EXDOSE[i]/myout$EXDOSE[i-1])-1)*100),1) ,3) #computing the increment to the next dose for the 0DLT
}
require(plyr)
myoutt <- ddply(myout,.(x),function(x){
x[1,]
})
return (myoutt)
}
#############################!!!!!!
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.
Date r params$d
knitr::opts_chunk$set(echo = TRUE,cache =FALSE)
# Multiple plot function # # ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) # - cols: Number of columns in layout # - layout: A matrix specifying the layout. If present, 'cols' is ignored. # # If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), # then plot 1 will go in the upper left, 2 will go in the upper right, and # 3 will go all the way across the bottom. # multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { library(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } }
### library(vcd) library(MASS) library("eha") library("msm") library("crmPack") source("safetywindow.r") source("mDA-CRM-model.r") source("mDA-CRM-design.r") library (xlsx) ##;
############ rules agreed with team ## Target toxicity interval: 20% – 35% DLT rate. ## Starting dose and dose range: 1.6 mg (to be confirmed) up to 400 mg (3 weekly flat dosing) ## Cohort size : at least 3 patients per cohort ## Maximum dose increments ## 100% increase prior to 1st DLT ## 50% increase after 1st DLT ## <25% probability to overdose (to be above targeted toxicity interval) ## Stopping criteria ## Chance (or likelihood) of being in the targeted toxicity interval >40% ## 6 subjects or at least 2 cohorts having dose within 20% dose of the MTD ## Clinical judgment: no benefit with additional cohorts ## A maximum of 60 patients is reached ## DLT window and safety windows ## DLT window 28 days from first administration. ## 1st patient safety window 1 week before enrolling the remaining patients, ## 1 day interval between 2nd and 3rd patient. ## for 2nd and 3rd patient at least 2 weeks safety window before starting the next cohort. ## Intra-patient escalation: allowed only if part of a dosing regimen explored, otherwise not planned ## from ceatcb safe dose it has been applied a 20 fold safety factor for determining the starting dose for cam5 ## starting dose for single patient cohort part is 65 microg ## doses of cea tcb are five fold higher than estimated dose of ceacam5 tcb so need to adjust to ## reflect the difference in potency myTmax <-28 #Tmax is the lenght of the DLT period in this project (CAM5) mynpiece <- c(myTmax/7) myStartDose <- 1.6 ### mydoseGrid=c(seq(from=0.5, to=10, by=0.1), seq(from=11, to=20, by=1),25,27, seq(from=30, to=100, by=5),120, seq(from=150, to=400, by=50)) set.seed(23)
################## Prior data from CEATCB ######################## ## Monotherapy study, N = 80 patients received at least one dose of IV RO6958688 QW up to 600 mg in ## escalation part of the trial (MTD for C1/C2 QW defined as 400 mg). Safety profile dominated by ## IRRs and Tumor Flare. ## (a dose grid up to 80 mg with 1 mg spacing, and from 85 to 1000 mg with 5 mg spacing,will be used). ## Five DLTs in five patients, at doses over 40 mg occurring shortly after the first administration ## 2 DLTs @ 600 mg (N = 2 pts), Respiratory failure Gr5, Colitis Gr4 (discontinued) ## 1 DLT @ 300 mg (N = 8 pts), Diarrhea Gr3 (resolved, C1D5-C1D18) ## 1 DLT @ 60 mg (N = 7 pts), Hypoxia Gr3 (resolved, C1D2-C1D10) ## 1 DLT @ 40mg (N = 12 pts), Dyspnea Gr3 (resolved, C1D2-C1D7) ## from protocol prior parameters for tcb model without pretreatment ## doseGrid 0.05 1 50 100 500 1000 ## mean 1.021, -0.682 ## cov 1.887,0.023,0.023,0.001 ## reference dose 500 #### reading in the tcb dataset dataTCB<- read.xlsx(file="dltdataTCB.xlsx",1) dataTCB$DLTDY[dataTCB$DLTDY>myTmax] <- myTmax ## note DLTDY cannot be bigger than the dlt window used later in the model dataTCB$DLTDY[dataTCB$DLTDY==0] <- 0.1 #looking at original set up for cea tcb mydoseGrid0 <- c(2.5, seq (3,85,1),seq(90,1000,5)) #from protocol myrefDose0<- c(500) model0 <- crmPack:::LogisticLogNormal(mean=c(1.021, -0.682), cov= matrix(c(1.887 , 0.023, 0.023, 0.001), nrow=2), refDose=myrefDose0) # keeping on the scale that is needed in this project #reading the data from ceatcb ; data0 <- crmPack:::Data(x=c(dataTCB$EXDOSE),y=c(dataTCB$DLT), doseGrid=mydoseGrid0) options <- crmPack:::McmcOptions(burnin=10000, step=2, samples=50000) set.seed(94) samples0 <- crmPack:::mcmc (data0,model0,options) plot(samples0,model0,data0) ### WARNING THE MTD IN HEREIS NOT AROUND 400.mtc was edfine also on the basis of other clinical consideration and the model prior was too informative ####
# Approach 1; # Based on the mixture prior to approximate a lognormal prior; ### getting mean and cov priors from mixture of informative and non informative components for ceatcb: myrefDose<- c(500/5) #scaling ref dose as well the doses from tcb dataset model_TCBS <- crmPack:::LogisticLogNormal(mean=c(1.021, -0.682), cov= matrix(c(1.887 , 0.023, 0.023, 0.001), nrow=2), refDose=myrefDose) # keeping on the scale that is needed in this project #reading the data from ceatcb scaling the doses; data_TCBS <- crmPack:::Data(x=c(dataTCB$EXDOSE/5),y=c(dataTCB$DLT), doseGrid=mydoseGrid) options <- crmPack:::McmcOptions(burnin=10000, step=2, samples=50000) set.seed(94) samples_TCBS <- crmPack:::mcmc (data_TCBS,model_TCBS,options) plot(samples_TCBS,model_TCBS,data_TCBS)
```{R inf_mod}
informative component
inmodel<- approximate(samples_TCBS, model_TCBS, data_TCBS, points = c(1,20,50,100,200,400), refDose = myrefDose, logNormal = TRUE, verbose = TRUE)
creating an empty dataset, provide only the dose grid:
emptydata_inm <- Data(doseGrid=mydoseGrid)
obtain prior samples with this Data object note inmodel is with posterior parameters !
priorsamples_inm <- mcmc(emptydata_inm, inmodel, options)
## then produce the plot plot(priorsamples_inm, inmodel, emptydata_inm) pi<- plot(priorsamples_inm, inmodel, emptydata_inm) #saving the plot as object to plot it later with multiplot
```r ## minimal inf component coarseGrid <- c(1,20,50,100,200,400) minInfModel <- MinimalInformative(dosegrid = coarseGrid, refDose=myrefDose, threshmin=0.1, threshmax=0.3, logNormal = TRUE) minInfModel$model@cov ## be carefyll that these two plots are one on the top of the other so check specification for output! matplot(x=coarseGrid, y=minInfModel$required, type="b", pch=19, col="blue", lty=1, xlab="dose", ylab="prior probability of DLT") matlines(x=coarseGrid, y=minInfModel$quantiles, type="b", pch=19, col="red", lty=1) legend("right", legend=c("quantiles", "approximation"), col=c("blue", "red"), lty=1, bty="n") noninmodel<-minInfModel$model emptydata <- Data(doseGrid=mydoseGrid) priorsamples_min <- mcmc(emptydata, noninmodel, options) ## then produce the plot plot(priorsamples_min, noninmodel, emptydata) pni<-plot(priorsamples_min, noninmodel, emptydata)#saving the plot as object to plot it later with multiplot
## mix informative and minimal informative models exploring different weights; components<- list(inf=list(mean=inmodel@mean,cov=inmodel@cov), noninf=list(mean=noninmodel@mean,cov=noninmodel@cov)) #for weight use the same order for the model as above, exploring different weights; mixmodel_1090 <- LogisticNormalFixedMixture(components=components, weights=c(10,90), refDose=myrefDose, logNormal = TRUE) mixmodel_2575 <- LogisticNormalFixedMixture(components=components, weights=c(25,75), refDose=myrefDose, logNormal = TRUE) mixmodel_4060 <- LogisticNormalFixedMixture(components=components, weights=c(40,60), refDose=myrefDose, logNormal = TRUE) mixmodel_7030 <- LogisticNormalFixedMixture(components=components, weights=c(70,30), refDose=myrefDose, logNormal = TRUE) emptydata <- Data(doseGrid=mydoseGrid) priorsamples_1090 <- mcmc(emptydata, mixmodel_1090, options) plot(priorsamples_1090, mixmodel_1090, emptydata) ## then produce the plot p1<- plot(priorsamples_1090, mixmodel_1090, emptydata) priorsamples_2575 <- mcmc(emptydata, mixmodel_2575, options) plot(priorsamples_2575, mixmodel_2575, emptydata) ## then produce the plot p2<- plot(priorsamples_2575, mixmodel_2575, emptydata) priorsamples_4060 <- mcmc(emptydata, mixmodel_4060, options) plot(priorsamples_4060, mixmodel_4060, emptydata) ## then produce the plot p3<- plot(priorsamples_4060, mixmodel_4060, emptydata) priorsamples_7030 <- mcmc(emptydata, mixmodel_7030, options) plot(priorsamples_7030, mixmodel_7030, emptydata) ## then produce the plot p4<- plot(priorsamples_7030, mixmodel_7030, emptydata) #more plots on one page multiplot(pi, pni, p1, cols=1) multiplot(p1, p2, p3, cols=1) multiplot(p1, p3, p4, cols=1) #overlaying plots p1innoin <- p1 + geom_line(data=pni$data, aes(x=x, y=y, group=group, linetype=Type), colour="blue" )+ geom_line(data=pi$data, aes(x=x, y=y, group=group, linetype=Type), colour="black") #make sure the caption is correct!! p1innoin+ labs (caption =c("In black the Inf. model, in blue the min. inf. mod. and in red the mixture with weight 10/90") ) p1234 <- p1 + geom_line(data=p2$data, aes(x=x, y=y, group=group, linetype=Type), colour="blue")+ geom_line(data=p3$data, aes(x=x, y=y, group=group, linetype=Type), colour="purple")+ geom_line(data=p4$data, aes(x=x, y=y, group=group, linetype=Type), colour="black") p1234 +labs (caption =c("Mixed models weights (inf/min.inf.): in red 10/90, in blue 25/75, in purple 40/60, in black 70/30") ) ## need to decide on a mixture prior for next steps and get the parameters ## get the 25/75 #overlaying plots p2innoin <- p2 + geom_line(data=pni$data, aes(x=x, y=y, group=group, linetype=Type), colour="blue" )+ geom_line(data=pi$data, aes(x=x, y=y, group=group, linetype=Type), colour="black") #make sure the caption is correct!! p2innoin+ labs (caption =c("In black the Inf. model, in blue the min. inf. mod. and in red the mixture with weight 25/75") )
posterior_hist <- approximate(priorsamples_2575, mixmodel_2575, emptydata, refDose=mixmodel_2575@refDose, logNormal = TRUE, seed = 12345, control = list(threshold.stop = 0.01, maxit = 50000, temperature = 50000, max.time = 120)) posterior_hist@mean posterior_hist@cov posterior_hist@refDose emptydata <- Data(doseGrid=mydoseGrid) priorsamples_post <- mcmc(emptydata, posterior_hist, options) plot_priorsam <- plot(priorsamples_post , posterior_hist, emptydata) multiplot(p2,plot_priorsam) #posterior_hist2 <- approximate(priorsamples_2575, mixmodel_2575, emptydata, #refDose=mixmodel_2575@refDose, logNormal = FALSE) #posterior_hist2@mean #posterior_hist2@cov #posterior_hist2@refDose #emptydata <- Data(doseGrid=mydoseGrid) #priorsamples_post2 <- mcmc(emptydata, posterior_hist2, options) #plot_priorsam2 <- plot(priorsamples_post2 , posterior_hist2, emptydata) # multiplot(p2,plot_priorsam,plot_priorsam2) ####Approach 2; ####Use quantiles to log normal function to estimate the mixture prior; # test <- Quantiles2LogisticNormal(dosegrid=c(1, 10, 30, 40, 80), # refDose=56, # lower= c(0.01,0.02,0.03,0.05,0.12), # median= c(0.02,0.08,0.3,0.45,0.77), #median probability of DLT at each dose in dosegrid is enough to put 5 to 6 doses, ## these are derived from the posterior logistic dose # upper= c(0.10,0.6,0.92,0.97,0.99), # logNormal=TRUE)
######################### modelling both dose and time ############################### ## reading in data including time of DLT MAKE SURE there are no 0 time!! data <- DataDA(x=c(dataTCB$EXDOSE)/5, ### patient dose level doses of cea tcb are five fold higher than estimated dose of ceacam5 tcb so need to adjust y=c(dataTCB$DLT), ### patient DLT 1=yes 0 =no doseGrid=mydoseGrid, u=c(dataTCB$DLTDY), ###u= DLT free survival values time to dlt or time to censoring if no dlt. ## NOTE u entries must not be larger than Tmax and ### dlt occuring outside the considered dlt window should be excluded. t0=c(dataTCB$dayd), ### time from recruitment of first patient to recruitment of each of the remaining patients Tmax=myTmax, npiece=mynpiece) #2) Structure of the model class #mean=c(-0.85,1) #cov=matrix(c(1,-0.5,-0.5,1),nrow=2) #refDose=56 #non-informative prior for lambdas (number of pieces=npiece_) #l=as.numeric(t(apply(as.matrix(c(1:npiece_),1,npiece_),2,lambda_prior))) ##need to fill in npiece_=mynpiece Tmax_=myTmax ## note dose grid in empty data need to be same as above emptydata_mda <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax,npiece=mynpiece) ### do not change the lambda prior is based on paper (is a non non informative lambda) lambda_prior<-function(k){ npiece_/(Tmax_*(npiece_-k+0.5)) } ### for lambda if want to use informative then may be better to fit from exponential otherwise may be the posterior of lambda is not stable enough ### coud try both with informative and with noninf lambda and see the model performance in the simulations ##try non informative first mymda_model<-DALogisticLogNormal(mean=c(posterior_hist@mean[1],posterior_hist@mean[2]), ##priors parameters use posterior derived from tcb mixture cov=matrix(c(posterior_hist@cov[1],posterior_hist@cov[2],posterior_hist@cov[3],posterior_hist@cov[4]),nrow=2), refDose=posterior_hist@refDose, l=as.numeric(t(apply(as.matrix(c(1:npiece_),1,npiece_),2,lambda_prior))), ##lambda here not informative and assuming higher prob of event later in the dlt time C_par=2) ## do not change C_par =2 it control the credible intervals around the H rate. ### ##priors parameters use posterior derived from tcb mixture #TRY DIFFRENT LAMBDA TO SEE THE EFFECT mymda_model2<-DALogisticLogNormal(mean=c(posterior_hist@mean[1],posterior_hist@mean[2]), cov=matrix(c(posterior_hist@cov[1],posterior_hist@cov[2],posterior_hist@cov[3],posterior_hist@cov[4]),nrow=2), refDose=posterior_hist@refDose, #l=as.numeric(t(apply(as.matrix(c(1:npiece_),1,npiece_),2,lambda_prior))), l=rev(c(0.04081633, 0.05714286, 0.09523810, 0.28571429)), ##lambda here assuming higher probability at the beginning of the dlt period C_par=2) ## do not change C_par =2 it control the credible intervals around the H rate.
options <- crmPack:::McmcOptions(burnin=10000, step=2, samples=50000) #3) Obtain the posterior using existing data for example #set.seed(94) #samples <- mcmc (data,mymda_model,options) #4) use ggmcmc to diagnose library(ggmcmc) #alpha0samples=get(samples,"alpha0") #print(ggs_traceplot(alpha0samples)) #print(ggs_autocorrelation(alpha0samples)) #5) plot the model fit #plot(samples, mymda_model,data,hazard=TRUE) # plot(samples, model,data,hazard=FALSE)#option FALSE does not work!
## PLOTTING PRIOR mDA emptydata <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax, npiece=mynpiece) Priorsamples <- mcmc (emptydata,mymda_model,options) alpha0samples=get(Priorsamples,"alpha0") print(ggs_traceplot(alpha0samples)) print(ggs_autocorrelation(alpha0samples)) #plot the model fit plot(Priorsamples, mymda_model,emptydata,hazard=TRUE) plotPrmda <- plot(Priorsamples, mymda_model,emptydata,hazard=TRUE) emptydata <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax, npiece=mynpiece) Priorsamples2 <- mcmc (emptydata,mymda_model2,options) alpha0samples=get(Priorsamples2,"alpha0") print(ggs_traceplot(alpha0samples)) print(ggs_autocorrelation(alpha0samples)) #plot the model fit plot(Priorsamples2, mymda_model2,emptydata,hazard=TRUE) plotPrmda2 <- plot(Priorsamples2, mymda_model2,emptydata,hazard=TRUE)
#6) Escalation rules ##need to fill in (use the same rule in the section 8 of "using the package crmPack: introductory examples") myIncrements <- IncrementsRelativeDLT(DLTintervals=c(0, 1), increments=c(1, 0.5)) myNextBest <- NextBestNCRM(target=c(0.2,0.35), overdose=c(0.35,1), maxOverdoseProb=0.25) mySize<- CohortSizeConst(size=3) myStopping1 <- StoppingTargetProb(target=c(0.2, 0.35), prob=0.4) myStopping2 <- StoppingMinPatients(nPatients=60) #should this be 60 instead?? myStopping3 <- StoppingPatientsNearDose(nPatients=6,percentage=20) myStopping <- (myStopping1 | myStopping2 | myStopping3) #7) recommended dose for the next cohort given the data ##### but data here are the tcb data..is this what is intended??is just an example of how it will behave? #nextMaxDose <- maxDose(myIncrements,data=data) # this depend on data!! #doseRecommendation <- nextBest(myNextBest, # doselimit=nextMaxDose, # samples=samples, ##need to run the mcm to get the posterior with current data # model=mymda_model, # data=data) #doseRecommendation$plot #doseRecommendation$value
##Example 2: run a simulation to evaluate design operating characters; #1) set up safety window and DADesign mysafetywindow=SafetyWindowConst(c(7,1),14,14) ####I could use SafetyWindowSize or SafetyWindowConst ### for each cohort I want to follow first patient for 7 days then all for at ### least 14 days #SafetyWindowConst(patientGap = c(7,3), # patientFollow = 7, this is the minimal required condition for all patienst in the cohort # patientFollowMin = 14) min is condition for at least one patient emptydata <- DataDA(doseGrid=mydoseGrid,Tmax=myTmax, npiece=mynpiece) design <- DADesign(model=mymda_model, increments=myIncrements, nextBest=myNextBest, stopping=myStopping, cohortSize=mySize, data=emptydata, safetyWindow=mysafetywindow, startingDose=c(myStartDose)) trial_bh <- examine(design,options=options) ## need to save output of trial_bh and run it overnight is slow!! ## some notes: ## DLTearly_1=1 indicates the scenario that the patients with longest follow up have DLTs, DLTearly_1=2 indicates The patients with shortest follow up have DLTs #question: it seems that at a certain point trial stop si always = true even in absence of dl and even if the next dose is higher than the current one why? #answer: In order to review more scenarios, the examine function stops only when newDose <= thisDose or the maximum dose in the dose grid reached. if stop=TRUE that means the stopping rule is fulfilled (this may due to a relatively informative prior). A quick way to check could be: write down the 0 DLT data (till 40mg for example) and run with the mcmc function. #Also why does it get to up to 6 dlt? are these the sum of the dlt on the current dose and the previous at teh end of the dlt period? #Yes, sum of the dlt on the current dose and the previous dose who are under the DLT follow up period, i.e. the nfollow I mentioned in the examine code. # also why even at the initial dose with 3 dlt keep going up? # This may due to a strong study prior. A way may confirm the reason is to run a nCRM model with the same alpha parameters. print(proc.time() - t1) #saving the file (is saving to the working directory) save(trial_bh,file="trial_bh.RData") trial_0 <- trial_bh[which (trial_bh$DLTs==0),] #setting a data set with 0 dlt assuming time of dosing between cohorts is as in safety window # need a data set with the following variable # $EXDOSE # $DLT # $DLTDY # $dayd ###################TO DO !!!!#####
#simulation with true scenarios #2)set up truth curves myTruth<-function(dose){ model@prob(dose,alpha0=2,alpha1=3) } curve(myTruth(x), from=0, to=c(max(mydoseGrid)), ylim=c(0, 1)) myTruth1<-function(dose){ model@prob(dose,alpha0=0.05,alpha1=2) } curve(myTruth1(x), from=0, to=c(max(mydoseGrid)), ylim=c(0, 1)) ### onset <- 15 ## median onset of dlt time feed in function below ## need to double check with the value of the exponential sample 1000 and see where is the mean mytrueTmax<- 28 #Truemax is used to generate the DLT data in the simulation and cannot be lower than onset, so any dlt is contained within 0 and mytrueTmax also feed in the function below. it could be set also bigger than the tmax used in the actual model to see what would happen if the dlt window used is in reality too short. # if the safety window is set as long as the dlt window then the simulation will evaluate a crm EWOC without mDA so then if the dlt window is shorter than truetmax it is possible to evaluate the advantage of da over other approache so the advantage of using in this case a longer dlt window but with short safety window exp_cond.cdf<-function(x){ 1-(pexp(x,1/onset,lower.tail=FALSE)-pexp(mytrueTmax,1/onset,lower.tail=FALSE))/pexp(mytrueTmax,1/onset) } #3) set up simulation settings mySims <- simulate(design, args=NULL, truthTox=myTruth, truthSurv=exp_cond.cdf,#piece_exp_cond.cdf, trueTmax=mytrueTmax, ### Truemax is used to generate the DLT data in the simulation. You can image if DLT is generated from #day 1 to day 28 in the simulation, but tmax=14 days, the design will not capture the full DLT information and it #leads to underestimation of the MTD. nsim=10, seed=819, mcmcOptions=options, firstSeparate=TRUE, deescalate=FALSE, ## it refer to patients already enrolled so with false if a dlt occur in a previous dose or current dose level the patients stay at the level they are. ## if set to TRUE if they deescalate they will be censored at the time of the occurrence of the dlt which cause the decision to de-escalate. parallel=FALSE) ## NOTE THE OPTION TRUE DOES NOT WORK YET FOR THIS mDA CRM FUNCTION!!!! # system.time(simulate(design, # args=NULL, # truthTox=myTruth, # truthSurv=exp_cond.cdf,#piece_exp_cond.cdf, # trueTmax=80, # nsim=500, # seed=819, # mcmcOptions=options, # firstSeparate=TRUE, # deescalate=FALSE, # parallel=FALSE)) #4) interpreate simulation result #use a similar way as section 9.2 in the "using the package crmPack: introductory examples" document a=summary(mySims,truth=myTruth) plot(mySims) mySims@stopReasons[[2]] savePlot <- function(myPlot,name) { png(filename = paste(Sys.Date(),"C:/Users/liaoz4/Documents/R/simulation_result/",name,".png",sep=""), width = 480, height = 480) print(myPlot) dev.off()}
###### this part is work on progress!! myOut<- function (model=mymda_model, Tmax=myTmax, npiece=mynpiece, doseGrid=mydoseGrid){ myout <- data.frame (i=c(1:60),EXDOSE=c(1.6,1.6,1.6,rep(NA,57)),DLT=rep(NA,60),DLTDY=c(rep(NA,60)) ,increment=rep(NA,60), x0DLT=c(1.6,1.6,1.6,rep(NA,57))) DLTDY=c(22,15,14) for (i in seq(4,58,3) ) rowSums(dat[,c("b", "c")], na.rm=TRUE) { #compute time from 1st patient in chortind <- (i-1)/3 timeind <- c(rep (DLTDY,each = 1, times = chortind)) myout$DLTDY[1:(i-1)][is.na(myout$DLTDY[1:(i-1)])] <- 0 myout$DLTDY[1:(i-1)] <- myout$DLTDY[1:(i-1)] + timeind[1:(i-1)] myout$DLTDY[myout$DLTDY>Tmax] <- Tmax #entries must not be larger than Tmax #computing the dose for the three patient cohort data <- DataDA(x=c(myout$EXDOSE), y=c(myout$DLT), ### patient DLT 1=yes 0 =no doseGrid=doseGrid, u=c(myout$DLTDY), ###u= DLT free survival values time to dlt or time to censoring if no dlt. ## NOTE u entries must not be larger than Tmax and ### dlt occuring outside the considered dlt window should be excluded. t0=c(myout$dayd), ### time from recruitment of first patient to recruitment of each of the remaining patients Tmax=Tmax, npiece=npiece) nextBest(myNextBest, doselimit=nextMaxDose, samples=samples, ##need to run the mcm to get the posterior with current data model=model, data=data) #myout$EXDOSE[i:(i+2)]<- rep( nextDose (myout$EXDOSE[1:i-1], myout$DLT[1:i-1], model),3) # to have the table in the column format myout$x0DLT[(i-3):(i-1)]<-myout$EXDOSE[i:(i+2)] myout$DLT[1:(i-1)] <- rep (0,c(i-1)) #reset dlt to 0 myout$increment[(i-3):(i-1)] <-rep( round((((myout$EXDOSE[i]/myout$EXDOSE[i-1])-1)*100),1) ,3) #computing the increment to the next dose for the 0DLT } require(plyr) myoutt <- ddply(myout,.(x),function(x){ x[1,] }) return (myoutt) } #############################!!!!!!
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.