source("/Users/pm/s/first.s")
dir <- "/Users/pm/book/ctrial/pm/3/"
setwd(dir)

##################################################################
# implements
# Ji, Li and Bekele, "Dose finding [..] based on toxicity intervals"
##################################################################

## doses
doses <- 1:8
ndoses <- length(doses) 
scenario1 <-  c(5,25,50,60,70,89,90,95)/100
xi <- 0.95; K1 <- 1; K2 <- 1.5    # design parameters
a0 <- 0.005; b0 <- 0.005          # prior pars
pT <- 0.20                        # target toxicity


p0 <- scenario1; nmax <- 30
sim.trial <- function(nmax=30,p0=scenario1)
  { ## simulate trial under scenario p0
    d <- 1      # initial dose
    y <- matrix(0,nrow=2,ncol=length(doses)) # data
                ## rows: n successes (row 1) and failures (2)
                ## cols: doses
    traj <- NULL ## record simuluation history
    k <- 3      # cohort size
    
    repeat{
      yi <- rbinom(1,k,p0[d])          # simulate response
      y[1,d] <- y[1,d]+yi              # update data
      y[2,d] <- y[2,d]+(k-yi)
      n <- sum(y)
      post <- posterior(y)             # update posterior
      traj <- rbind(traj,c(d,post$m))  # update history 
      d <- dose.assign(d,y,post)       # dose assignment for next cohort
        ## returns 0 to indicate stopping if lowest dose is too toxic
      if ( (n>nmax) | (d==0) )
        break
    }
    colnames(traj) <- c("DOSE",paste("M",doses,sep=""))
    ds <- dose.dlt(y)               # final DLT estimation
    mx.dose <- ds$mx.dose
    return(list(y=y, dstar=ds$dstar, m=ds$m, traj=traj[,1:(1+mx.dose)]))
  }

posterior <- function(y)
  { ## posterior updates
    a1 <- y[1,]+a0; b1 <- y[2,]+b0        # posterior parameters, p[d]~Be(a1[d],b1[d])
    m <- a1/(a1+b1)                       # moments of p[d]
    s <- sqrt( a1*b1/(a1+b1+1) )/(a1+b1)  
    return(list(a1=a1,b1=b1,m=m,s=s))
  }

dose.assign <- function(d,y,post)
  { ## find optimal dose for next cohort, starting from current dose d
    ## returns d=0 if dose 1 is too toxic
    pExcess <- 1-pbeta(pT, post$a1,post$b1)          # p(pi>pT | y) (needed later)
    if ( (d==1) & (pExcess[1]>xi))
      dnew <- 0
    else{
      s <- post$s[d]; a1 <- post$a1[d]; b1 <- post$b1[d]
      ## probabilities of the three interval partition
      qD <- 1-pbeta(pT+K1*s, a1,b1) # p(pi > pT+K1*s | y)
      qE <-   pbeta(pT-K2*s, a1,b1) # p(pi < pT-K2*s | y)
      qS <- 1-(qD+qE)                        # p(pT-K2*s < pi < pT+K2*s | y)
      ## record utilities of possible actions (D,S,E) using u= -1 for impossible action
      u <- c(qD,qS,qE)              # initialize utilities
      if (d<ndoses){
        if (pExcess[d+1]>xi)        # no escalation if toxcitity is too high  
          u[3] <- -1
      }
      if (d==1) u[1] <- -1
      dnew <- d+which.max(u)-2 
      if ( (dnew<1) | (dnew>ndoses) )
        cat("\n *** error: impossible dose?\n")
    }
    return(dnew)
  }
    
dose.dlt <- function(y)
  { ## report DLT
    a1 <- y[1,]+a0; b1 <- y[2,]+b0        # posterior parameters, p[d]~Be(a1[d],b1[d])
    m <- a1/(a1+b1)                       # moments of p[d]
    ## restrict search to below max assigned dose
    mx.dose <- last(which(apply(y,2,sum)>0))
    m1 <- monotone(m[1:mx.dose])
    cat("Selecting max dose from 1...",mx.dose,"\n")
    dstar <- which.min(abs(m1-pT))
    return(list(dstar=dstar,m=m1,mx.dose=mx.dose))
  }

monotone <- function(m)
  { # pool adjcacent violators
    p <- length(m)
    idx <- 1:p # equal indices indicate pooling
      ## initialize w/o any pooling
    m1 <- m
    while(any(diff(m1)<0)){
      i <- which(m1[-p]>m1[-1])[1] # find 1st violator
      i1 <- is.element(idx, idx[i:(i+1)]) # pool the two adjacent guys
      m1[i1] <- mean(m[i1])
      idx[i1] <- min(idx[i1])
    }
    return(m1)
  }


example <- function()
  { ## execute these lines for example
    st <- sim.trial(p0=scenario1)

    M <- nrow(st$traj)
    ps("324-traj")
    matplot(1:M,st$traj[,-1],type="l",xlab="WEEK",
            ylab="E(TOX | dose,  Y)",bty="l")
    nd <- ncol(st$traj)-1
    legend(M-4,0.6,legend=c("DOSE",1:nd),lty=0:nd,col=0:nd,bty="n")
    devoff()

    ps("324-d")
    plot(1:M,st$traj[,1],type="p",pch=17,xlab="WEEK",ylab="DOSE",
         bty="l")
    abline(h=st$dstar,col=2,lty=3)
    text(0.71,st$dstar-0.1 ,adj= 0,"D*",col=2)
    devoff()
  }