## Dec 9 2006, Yuan Ji
## This is the R code for perform phase I dose-finding method in Ji, Li and Bekele (2007), which uses independent beta priors. The results are written to the file "simresult.txt". The end users have to provide the simulation parameters "simN", "pT", "D", and "sampsize", which are explained below. Also, the scenarios that determine the true dose probabilities must be provided. This program uses the 10 scenarios in Ji, Li and Bekele (2007)

##library(SAGx)
system("rm -f simresult.txt")

## pava is the pool adjacent violator algorithm to perform isotonic transformation for the posterior means later

pava <- function (x, wt = rep(1, length(x))) 
{
    n <- length(x)
    if (n <= 1) 
        return(x)
    if (any(is.na(x)) || any(is.na(wt))) {
        stop("Missing values in 'x' or 'wt' not allowed")
    }
    lvlsets <- (1:n)
    repeat {
        viol <- (as.vector(diff(x)) < 0)
        if (!(any(viol))) 
            break
        i <- min((1:(n - 1))[viol])
        lvl1 <- lvlsets[i]
        lvl2 <- lvlsets[i + 1]
        ilvl <- (lvlsets == lvl1 | lvlsets == lvl2)
        x[ilvl] <- sum(x[ilvl] * wt[ilvl])/sum(wt[ilvl])
        lvlsets[ilvl] <- lvl1
    }
    x
}

## Need to specify these simulation parameters before run the simulations
simN<-1000  ## 1000 simulations

pT<-.2; D<-4; sampsize <- 18; csize <-2; startdose <- 2
## target tox probility is pT, with D doses, cohorts size is csize, and max sample size is sampsize


set.seed(2)

betavar<-function(a,b){
  resp <- a*b/((a+b)^2*(a+b+1))
  return(resp)
}

## These parameters are recommended to be fixed for non-advanced users.
a<-.005; b<-.005;  xi<-.95
k1<-1; k2<-1.5; pow <- 1 ## pow is the alpha in the model



## Set up scenarios (the program include 10 scenarios that are in Ji, Li and Bekele.
for(sc in c(1:10)){
  if (sc==1){
    p <- c(0.05, 0.25, .5, .6 ,.7, .8, .9, .95)
  }
  if (sc==2){
    p <- c(0.01, 0.05, .25, .5 ,.6, .7, .8, .9)
  }
  if (sc==3){
    p <- c(.01, .02, .05, .25, .5, .6, .7, .8)
  }
  if (sc==4){
    p <- c(.01, .02, .03, .05, .25, .5, .6, .7)
  }
  if (sc==5){
    p <- c(.01, .02, .03, .04, .05, .25, .5, .6)
  }
  if (sc==6){
    p <- c(.01, .02, .03, .04, .05, .06, .25, .5)
  }
  if (sc==7){
    p <- c(.01, .05, .5, .6, .7, .8, .9, .95)
  }
  if (sc==8){
    p <- c(.4, .5, .6, .7, .8, .9, .95, .99)
  }
  if (sc==9){
    p <- c(.15, .25, .35, .45, .55, .65, .75, .85)
  }
  if (sc==10){
    p <- c(.05, .15, .25, .35, .45, .55, .65, .75)
  }
 

datan<-matrix(rep(0,simN*D),ncol=D)
datax<-matrix(rep(0,simN*D),ncol=D)

rez<-rep(0,simN)

for(sim in 1:simN){
  ##print(sim)
  x<-rep(0,D); n<-rep(0,D)
  pa<-rep(0, D); pb<-rep(0,D)
  q <- rep(0,3)
  d<-startdose; st<-0; nodose<-0; maxdose<-1; toxdose<-D+1; seldose<-0
  
  
  ## Start simulation; Generate data
  while(st==0){
    maxdose<-max(maxdose, d)
    xx <- 0
    for(i in 1:csize){
      ttt <- runif(1)
      if(ttt < p[d]) xx <- xx + 1
    }
    ##xx<-rbinom(1, 3, p[d])
    x[d] <- x[d] + xx; n[d] <- n[d] + csize
	a<-.005; b<-.005
	
    pa[d]<-x[d]+a; pb[d]<-n[d]-x[d]+b
    pa[d+1]<-x[d+1]+a; pb[d+1]<-n[d+1]-x[d+1]+b
    
    ##Compute the indicator T_{i+1}
    if(d<D){
      temp<-1-pbeta(pT, pa[d+1], pb[d+1])
      if(temp>xi) {tt<-1; toxdose<-d+1} else {tt<-0}
    }
	
    ##Compute the three posterior probabilities
    postsd <- sqrt(betavar(pa[d], pb[d]))
    q[1] <- 1-pbeta(k1*postsd^(pow)+pT, pa[d], pb[d])
    q[2] <- pbeta(k1*postsd^(pow)+pT, pa[d], pb[d]) - pbeta(pT-k2*postsd^(pow), pa[d], pb[d])
    q[3] <- (1-q[1]-q[2])*(1-tt)
    
    ## implement the dose-assignment rules
    if(d==1){
      ## if the first dose is too toxic, the trial will be terminated
      if((1-pbeta(pT, pa[d], pb[d]))>xi){st=1; nodose<-1;} 
      else{
        if((q[2]>q[1])&&(q[2]>q[3])){d<-d}
        if((q[3]>q[1])&&(q[3]>q[2])){d<-d+1}
      }
    }
    else{
      if(d==D){
        if((q[1]>q[2])&&(q[1]>q[3])){d<-d-1}
        if((q[2]>q[1])&&(q[2]>q[3])){d<-d}
      }
      else{
        if((d>1)&&(d<D)){
          if((q[1]>q[2])&&(q[1]>q[3])){d<-d-1}
          if((q[2]>q[1])&&(q[2]>q[3])){d<-d}
          if((q[3]>q[1])&&(q[3]>q[2])){d<-d+1}
        }
      }
    }
    total<-sum(n)
    if(total >= sampsize){st<-1}
  }
  
  ## compute the posterior mean from the beta distribution
  if(nodose==0){
    tdose<-min(maxdose, toxdose-1)
    pp<-rep(-100,tdose)
    pp.var<-rep(0, tdose)
    
    for(i in 1:tdose){
      if(pa[i]!=0){pp[i] <- pa[i]/(pa[i]+pb[i]); pp.var[i] <- betavar(pa[i], pb[i])}
    }
  
    pp<-pava(pp, wt=1/pp.var)  ## perform the isotonic transformation using PAVA with weights being posterior variances
	  
    for(i in 2:tdose){
	 #pp[i] <- max(pp[i], pp[i-1]+.00001) 
	 pp[i] <- pp[i] + i*1E-10 ## by adding an increasingly small number to tox prob at higher doses, it will break the ties and make the lower dose level the MTD if the ties were larger than pT or make the higher dose level the MTD if the ties are smaller than pT
	}
    seldose<-sort(abs(pp-pT), index.return=T)$ix[1]
    ##seldose is the final MTD that is selected
  }
  rez[sim] <- seldose;
  for(i in 1:D){
    datan[sim,i] <- n[i]
    datax[sim,i] <- x[i]
  }
}
##rez[is.na(rez)]<-0
aaa<-rep(0,D)

for(i in 1:D){
  aaa[i] <- sum(rez==i)/simN
}
write(sc, "simresult.txt", append=T) ## Scenario number
write(p, "simresult.txt", ncolumns = 8, append=T) ## the true tox probabilities
write(aaa, ncolumns = 8, "simresult.txt", append=T) ## aaa is the vector containing the percentages of selection
write(apply(datan,2,mean), ncolumns = 8, "simresult.txt", append=T) ## prints out the average number of patients at each dose
write(sum(apply(datan,2,mean)),  "simresult.txt", append=T) ## prints out the average total sample size
write(sum(datax)/sum(datan), "simresult.txt", append=T)  ## prints ot the overall toxicity
write("**********************", "simresult.txt", append=T)
}