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

##################################################################
# implements
# Thall, Simon & Estey (1995, Stat in Med)
##################################################################

## doses
K <- 4           ## number of elementary events (sec 2.1. of Thall & al)
thS <- c(2.037,6.111,30.555,2.037)  ## pars for Dir prior for event prob's
                 ## table 1, Thall & al
thE <- thS/(sum(thS))*K   ## same for experimental therapy
CR <- 1:2        ## CR = A1 + A2
TOX <- c(2,4)    ##      A2 + A4
etaS <- 0.2     ## efficacy prob under standard of care
etaST <- 0.2     ## tox prob un
deltaCR <- 0.20  ## min improvement in eff
deltaTOX <- 0.05  ## slip factor for design for tox monitoring
pLC <-  0.02     ## lower threshold for improvement on efficacy resp
pUC <- 1.0       ## upper ... (not used in the example)
pUT <- 0.80
n.mx <- 75

## scenario 1:
make.p0 <- function(p12=0.20, p24=0.10)
  {# aux function to return 4 cell prob's assuming independence
    p2 <- p12*p24 # assume indep
    p1 <- p12-p2
    p4 <- p24-p2
    p3 <- 1-p1-p2-p4
    return(c(p1,p2,p3,p4))
  }
scenario1 <- make.p0(.2, .1)
## scenario 2:
scenario2 <- make.p0(.2,.2)
scenario3 <- make.p0(.2,.3)
scenario4 <- make.p0(.4,.1)
p0 <- scenario4

sim.trial <- function(n.mx=75,p0=scenario1)
  { ## simulate trial under scenario p0
    th1 <- thE # posterior parameters, initialize post=prior
    traj <- NULL
    y <-  NULL    ## patient responses and recruitment times (in months)
    month <- 0    ## calendar time -- need to deterine 3 month horizon
    n1 <- 0       ## number observed patients
    n <-  0       ## n recruited patients
    repeat{
      p12 <- post.probs(th1)  # post probs of events of interest
      cat("n=",n1,p12,"\n")
      traj <- rbind(traj,p12)
      if (p12[1] <= pLC){
        cat("stopping trial for lack of efficacy.\n")
        break
      } else if (p12[2] >= pUT) {
        cat("stopping trial for excess toxicity.\n")
        break
      } else if (n1>=n.mx){
        cat("stopping trial for n max.\n")
        break
      }
      ## recruit one more cohort of 4 patients per month
      ## and update post parameters
      if (n < n.mx){
      ## NOTE: continue trial until all patients have 3 months
      ##       i.e., last 2 months don't accrue anymore
        yi <- sample(1:4,4,prob=p0,replace=T)
        y <- rbind(y,cbind(yi,month))   # update data
        n <- n+4                        # update n recruited patients
      }
      month <-  month +1
      i0 <- which( y[,2]+3 <= month )   # response observed?
      n1 <- length(i0)                  # n recorded patients
      for (i in i0){
        yi <- y[i,1]
        th1[yi] <- th1[yi]+1            # update posterior parameters
      }
    }

    return(list(traj=traj,y=y))
  }

post.probs <- function(th1)
  { ## compute posterior probs of events of interest by MC simulation
    ## th1 = par of Dir posterior for experimental therapy
    ## thS = par of Dir prior for standard theray (doesn't change)
      M <- 100
      a1 <- sum( th1[CR] )    # experimental therapy
      b1 <- sum( th1[-CR] )
      a1S <- sum( thS[CR] )    # standard therapy
      b1S <- sum( thS[-CR] )
      a2 <- sum(th1[TOX])
      b2 <- sum(th1[-TOX])
      a2S <- sum(thS[TOX])
      b2S <- sum(thS[-TOX])
      eta1S <- rbeta(M,a1S,b1S)
      eta2S <- rbeta(M,a2S,b2S)
      p1 <- mean( 1 - pbeta(eta1S + deltaCR,  a1, b1) )
      p2 <- mean( 1 - pbeta(eta2S + deltaTOX, a2, b2) )
      return(c(p1,p2))
    }

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

    M <- nrow(st$traj)
    ps("431-traj")
    matplot(1:M,st$traj,type="l",xlab="MONTH",
            ylab="Pr(eta[S]+delta < eta[E] | Y)",bty="l")
    abline(h=c(.02,.8),col=1:2,lty=3)
    legend(M,0.6,xjust=1,legend=c("EFF","TOX"),lty=1:2,col=1:2,bty="n")
    devoff()

    ps("431-d")
    plot(1:M,st$traj[,1],type="p",pch=17,xlab="WEEK",ylab="DOSE",
         bty="l")
    devoff()
  }