#################
# duration.plot #
#################

duration.plot <- function(n,T,P,m,tm,B=1000,Tmax=2*T,sample.paths=0) {
# n    is the target sample size
# T    is the target completion time
# P    is the prior certainty (0 <= P <= 1)
# m    is the sample observed to date
# tm   is the time to date
# B    is the number of simulated duration times
# Tmax is the upper limit on the time axis of the graph
#
# set P to zero for a non-informative prior
# set m and tm to zero if no data has been accumulated yet.
  
  k <- n*P+m
  V <- T*P+tm
  theta <- 1/rgamma(B,shape=k,rate=V)
  simulated.duration <- matrix(NA,nrow=n-m,ncol=B)
  
  for (i in 1:B) {
    wait <- rexp(n-m,1/theta[i])
    simulated.duration[,i] <- tm+cumsum(wait)
  }
  
  lcl <- apply(simulated.duration,1,quantile,probs=0.025)
  mid <- apply(simulated.duration,1,quantile,probs=0.500)
  ucl <- apply(simulated.duration,1,quantile,probs=0.975)
  
  layout(matrix(c(1,2,2,2)))
  par(mar=c(0.1,4.1,0.1,0.1))
  duration.hist <- cut(simulated.duration[n-m,],seq(0,Tmax,length=40))
  barplot(table(duration.hist),horiz=FALSE,axes=F,xlab=" ",ylab=" ",space=0,col="white",names.arg=rep(" ",39))
  
  par(mar=c(4.1,4.1,0.1,0.1))
  plot(c(0,Tmax),c(0,n),xlab="Time",ylab="Number of patients",type="n")
  polygon(c(lcl,rev(ucl)),c((m+1):n,n:(m+1)),density=-1,col="gray",border=NA)
  lines(mid,(m+1):n,col="white")
  segments(0,0,tm,m)
  
  while (sample.paths>0) {
    lines(simulated.duration[,sample.paths],(m+1):n)
    sample.paths <- sample.paths-1
  }  
  
  pctiles <- matrix(NA,nrow=2,ncol=3)
  dimnames(pctiles) <- list(c("Waiting time","Trial duration"),c("2.5%","50%","97.5%"))
  pctiles[1,] <- 1/qgamma(c(0.975,0.5,0.025),shape=k,rate=V)
  pctiles[2,1] <- lcl[n-m]
  pctiles[2,2] <- mid[n-m]
  pctiles[2,3] <- ucl[n-m]
  list(pct=pctiles,sd=simulated.duration)
}
p0 <- duration.plot(n=350,T=3,P=0.5,m= 0,tm=      0,Tmax=10)
p1 <- duration.plot(n=350,T=3,P=0,  m=41,tm=239/365,Tmax=10)
p2 <- duration.plot(n=350,T=3,P=0.5,m=41,tm=239/365,Tmax=10)
p0
p1
p2

################
# accrual.plot #
################

accrual.plot <- function(n,T,P,m,tm,B=1000,nmax=2*n,sample.paths=0) {
# n    is the target sample size
# T    is the target completion time
# P    is the prior certainty (0 <= P <= 1)
# m    is the sample observed to date
# tm   is the time to date
# B    is the number of simulated accrual times
# nmax is the upper limit on the sample size axis of the graph
#
# set P to zero for a non-informative prior
# set m and tm to zero if no data has been accumulated yet.
  
  k <- n*P+m
  V <- T*P+tm
  theta <- 1/rgamma(B,shape=k,rate=V)
  simulated.duration <- matrix(NA,nrow=nmax-m,ncol=B)

  for (i in 1:B) {
    wait <- rexp(nmax-m,1/theta[i])
    simulated.duration[,i] <- tm+cumsum(wait)
  }

  tlist <- seq(tm,T,length=100)
  accrual.count <- function(x,t) {m+sum(x<=t)}
  simulated.accrual <- matrix(NA,nrow=100,ncol=B)

  for (i in 1:100) {
   time <- tlist[i]
   simulated.accrual[i,] <- apply(simulated.duration,2,accrual.count,t=time)
  }

  lcl <- apply(simulated.accrual,1,quantile,probs=0.025)
  mid <- apply(simulated.accrual,1,quantile,probs=0.500)
  ucl <- apply(simulated.accrual,1,quantile,probs=0.975)
  
  layout(matrix(c(2,2,2,1),nrow=1))
  par(mar=c(4.1,0.1,0.1,0.1))
  accrual.hist <- cut(simulated.accrual[100,],seq(0,nmax,length=40))
  barplot(table(accrual.hist),horiz=TRUE,axes=F,xlab=" ",ylab=" ",space=0,col="white",names.arg=rep(" ",39))

  par(mar=c(4.1,4.1,0.1,0.1))
  plot(c(0,T),c(0,nmax),xlab="Time",ylab="Number of patients",type="n")
  polygon(c(tlist,rev(tlist)),c(lcl,rev(ucl)),density=-1,col="gray",border=NA)
  lines(tlist,mid,col="white")
  segments(0,0,tm,m)

  while (sample.paths>0) {
    lines(tlist,simulated.accrual[,sample.paths])
    sample.paths <- sample.paths-1
  }  
  
  pctiles <- matrix(NA,nrow=2,ncol=3)
  dimnames(pctiles) <- list(c("Waiting time","Trial accrual"),c("2.5%","50%","97.5%"))
  pctiles[1,] <- 1/qgamma(c(0.975,0.5,0.025),shape=k,rate=V)
  pctiles[2,1] <- lcl[100]
  pctiles[2,2] <- mid[100]
  pctiles[2,3] <- ucl[100]
  list(pct=pctiles,sa=simulated.accrual)
}
d0 <- accrual.plot(n=350,T=3,P=0.5,m= 0,tm=      0,nmax=500,sample.paths=5)
d1 <- accrual.plot(n=350,T=3,P=0,  m=41,tm=239/365,nmax=500)
d2 <- accrual.plot(n=350,T=3,P=0.5,m=41,tm=239/365,nmax=500)
d0
d1
d2

##########################################################
# Example                                                #
#  P=0.5, k=175, V=1.5 years (547.5 days)                #
#  T=3 years (1095 days), n=350, T/n = 3.1 days          #
#  We observe m=41 patients, accrued at time tm=239 days #
#  tm/m = 5.9 days                                       #
##########################################################

p0 <- duration.plot(n=350,T=3,P=0.5,m= 0,tm=      0,Tmax= 10,sample.paths=2)
p1 <- duration.plot(n=350,T=3,P=0.0,m=41,tm=239/365,Tmax= 10,sample.paths=2)
p2 <- duration.plot(n=350,T=3,P=0.5,m=41,tm=239/365,Tmax= 10,sample.paths=2)
d0 <-  accrual.plot(n=350,T=3,P=0.5,m= 0,tm=      0,nmax=500,sample.paths=2)
d1 <-  accrual.plot(n=350,T=3,P=0.0,m=41,tm=239/365,nmax=500,sample.paths=2)
d2 <-  accrual.plot(n=350,T=3,P=0.5,m=41,tm=239/365,nmax=500,sample.paths=i)

###############
# end of file #
###############