In Roche/crmPack: Object-Oriented Implementation of CRM Designs

Date r params$d

knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
# nolint start

# 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

rv <- (approximate(samples_TCBS, model_TCBS, data_TCBS, points = c(1, 20, 50, 100, 200, 400), refDose = myrefDose, logNormal = TRUE, verbose = TRUE ))$model

rv$model

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"))

rv <- 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
  )
)
prior_hist < rv$model

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)

rv <- approximate(priorsamples_2575, mixmodel_2575, emptydata), #refDose=mixmodel_2575@refDose, logNormal = FALSE)
prior_hist2 <- rv$model

# 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(
  intervals = 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(
#   gap = c(7,3),
#   follow = 7, # this is the minimal required condition for all patienst in the cohort
#   follow_min = 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,
  cohort_size = 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)
}



############################# !!!!!!

# nolint end