R/boundcode_equalrand_jsm.R

In GBOP2: Generalized Bayesian Optimal Phase II Design (G-BOP2)

##library(Rcpp)
##library(RcppArmadillo)
##library(Rcpp11)
#library(MASS)
##library(RcppEigen)
##library(progress)
# sourceCpp("Calculation.cpp")
# sourceCpp(file="Calculation.cpp",cacheDir="cache")

maxresp_dual <- function(cutoff1,cutoff2,cutoff2_lambda1,cutoff2_lambda2,dprior,ntr,phi,delta0,delta1,contrast)
{ ##dprior: Dirichlet prior
  ##contrast: design vector: b with elements of 0 and 1
  ##ntr=n=nobs.seq[i]=nobs[i]:  sample sizes
  ##ntr total number of patients treated at current,0:ntr means number of repsonses from 0 to ntr.
  ncut_fut=NULL
  ncut_grad=NULL
  contrast = matrix(contrast,nrow=length(phi))

  ncut_fut = max(which(1-pbeta(phi+delta0,0:ntr+sum(dprior[which(contrast[1,]==1)]),ntr:0+sum(dprior)-sum(dprior[which(contrast[1,]==1)]))<(1-cutoff1)))-1
  if (1-cutoff2 != 0) {
    tmp = which(1-pbeta(phi+delta1 * (2-cutoff2),0:ntr+sum(dprior[which(contrast[1,]==1)]),ntr:0+sum(dprior)-sum(dprior[which(contrast[1,]==1)]))>1-cutoff2_lambda1)
    if (sum(tmp) == 0){
      ncut_grad = 999
    } else {
      ncut_grad = min(tmp)-1
    }
  } else {
    tmp = which(1-pbeta(phi,0:ntr+sum(dprior[which(contrast[1,]==1)]),ntr:0+sum(dprior)-sum(dprior[which(contrast[1,]==1)]))>cutoff2_lambda2)
    if (sum(tmp) == 0) {
      ncut_grad = 999
    } else {
      ncut_grad = min(tmp)-1
    }
  }

  results = c(ncut_fut, ncut_grad)
  return(results)
}

##Binary Efficacy, Ordinal Efficacy, Multiple Efficacy
# getboundary <- function(dprior,contrast,nobs,b,b2,b_grad1,b_grad2,pow,phi,delta0,delta1)
getboundary_dual <- function(dprior,contrast,nobs,b,b_grad1,b_grad2,pow1,pow2,pow3,phi,delta0,delta1)
{
  stopbound_fut = NULL
  stopbound_grad = NULL
  # nm = sort(union(nobs,nobsTX))
  nmax=max(nobs)
  for (n in nobs){
    # if (n == nmax){
    #   cutoffs = maxresp(cutoff1=b*(n/nmax)^pow1,cutoff2=(n/nmax)^(pow2),cutoff2_lambda=b_grad2,
    #                     dprior=dprior,contrast=contrast,ntr=n,phi=phi,delta0=delta0,delta1=delta1)
    # } else {
    #   ccutoffs = maxresp(cutoff1=b*(n/nmax)^pow1,cutoff2=(n/nmax)^(pow2),cutoff2_lambda=b_grad1,
    #                      dprior=dprior,contrast=contrast,ntr=n,phi=phi,delta0=delta0,delta1=delta1)
    # }
    # fff(b*pow(n/nmax,pow1),pow(n/nmax,pow2),b_grad1,b_grad2,dprior,n,phi,delta0,delta1,contrast);
    cutoffs = maxresp_dual(cutoff1=b*(n/nmax)^pow1,cutoff2=(n/nmax)^(pow2),
                      cutoff2_lambda1=b_grad1*(n/nmax)^(pow3),cutoff2_lambda2=b_grad2,
                     dprior=dprior,contrast=contrast,ntr=n,phi=phi,delta0=delta0,delta1=delta1)

    # message(cutoffs)

    stopbound_fut = cbind(stopbound_fut,cutoffs[1])
    if(!(n %in% nobs)) stopbound_fut[ncol(stopbound_fut)] = -Inf
    stopbound_grad = cbind(stopbound_grad,cutoffs[2])
    if(!(n %in% nobs)) stopbound_grad[ncol(stopbound_grad)] = -Inf
  }
  boundary <- rbind(nobs,stopbound_fut,stopbound_grad)
  boundary[which(boundary>999)]=999
  boundary[which(boundary< -999)]=-999

  return(boundary)
}

##bound of Binary Toxicity
getboundary.bitox <- function(dprior,contrast=as.matrix(1),nobs,b,pow,phi)
{
  stopbound = NULL
  nmax=max(nobs)
  for (n in nobs)
  {
    stopbound = cbind(stopbound,maxresp(cutoff=b*(n/nmax)^(pow/3),dprior=dprior,contrast=contrast,ntr=n,phi=phi))
  }
  boundary <- rbind(nobs,nobs-stopbound)
  boundary[which(boundary>999)]=999
  boundary[which(boundary< -999)]=-999
  ###patients who do not experience toxitity
  return(boundary)
}

##bound of Efficacy & Toxicity
#dimnames(bound.1) <- list(c("# patients treated","Stop if # CR <=","Stop if # CR/PR <="),NULL)
getboundary.tox <- function(dprior,contrast,nobs,nobsTX,b,pow,phi)
{
  stopbound1 = NULL
  stopbound2 = NULL
  nm = sort(union(nobs,nobsTX))
  nmax=max(nm)
  for (n in nm)
  {
    stopbound1 = cbind(stopbound1,maxresp(cutoff=b*(n/nmax)^pow,dprior=dprior,contrast=as.matrix(contrast[1,]),ntr=n,phi=phi[1]))
    if(!(n %in% nobs)) stopbound1[ncol(stopbound1)] = -Inf
    stopbound2 = cbind(stopbound2,maxresp(cutoff=b*(n/nmax)^(pow/3),dprior=dprior,contrast=as.matrix(contrast[2,]),ntr=n,phi=phi[2]))
    if(!(n %in% nobsTX)) stopbound2[ncol(stopbound2)] = -Inf

    # stopbound1 = cbind(stopbound1,maxresp(cutoff=b*(n/nmax)^pow,dprior=dprior,contrast=contrast,ntr=n,phi=phi))

  }
  # boundary <- rbind(nm,stopbound1[1,],nm-stopbound1[2,])

  boundary <- rbind(nm,stopbound1,nm-stopbound2)
  boundary[which(boundary>999)]=999
  boundary[which(boundary< -999)]=-999
  ##first row efficacy
  ## second row toxicity
  return(boundary)
}

getboundary.TTE <- function(gprior, phi, coeff, pow, n.interim, n.column=11)
{
  nmax = max(n.interim)
  boundary.full = c()
  #boundary.full = vector("list",length=length(n.interim))
  for (k in 1:length(n.interim))
  {
    boundary.temp = c()
    for (ss in 0:n.interim[k])
    {

      st = 0
      t.fu = 0
      while(st==0)
      {
        #if(pinvgamma(phi/log(2), shape = gprior[1]+ss, rate = gprior[2]+t.fu) <= 1 - coeff * (n.interim[k]/nmax)^pow) st=1

        if(pgamma(log(2)/phi, shape = gprior[1]+ss, rate = gprior[2]+t.fu,lower.tail = FALSE) <= 1 - coeff * (n.interim[k]/nmax)^pow) st=1
        else t.fu=t.fu+1
      }
      if(t.fu>0) {
        t.seq = seq(t.fu-1,t.fu,by=0.03)
        #t.bound = t.seq[max(which(pinvgamma(phi/log(2), shape = gprior[1]+ss, rate = gprior[2]+t.seq) > 1 - coeff * (n.interim[k]/nmax)^pow))]

        t.bound = t.seq[max(which(pgamma(log(2)/phi, shape = gprior[1]+ss, rate = gprior[2]+t.seq,lower.tail = FALSE) > 1 - coeff * (n.interim[k]/nmax)^pow))]
      }
      else t.bound=0

      boundary.temp = cbind(boundary.temp,c(n.interim[k],ss,t.bound))
    }
    boundary.full = cbind(boundary.full,boundary.temp)
  }
  return(boundary.full)
}


resp_list_r=function(seed, nsim,nobs,ptrue){

set.seed(seed)

resplist = list()
resplist_1=matrix(nrow = length(nobs), ncol = nsim)
resplist_2=matrix(nrow = length(nobs), ncol = nsim)
resplist_3=matrix(nrow = length(nobs), ncol = nsim)
resplist_4=matrix(nrow = length(nobs), ncol = nsim)
resplist[[1]] <- rmultinom(nsim,size=nobs[1],prob=ptrue)
resplist_1[1,] <- resplist[[1]][1,]
resplist_2[1,] <- resplist[[1]][2,]
resplist_3[1,] <- resplist[[1]][3,]
resplist_4[1,] <- resplist[[1]][4,]
for (j in 2:length(nobs))
{
  resplist[[j]] <- rmultinom(nsim,size=nobs[j]-nobs[j-1],prob=ptrue)
  resplist_1[j,]=resplist[[j]][1,]
  resplist_2[j,]=resplist[[j]][2,]
  resplist_3[j,]=resplist[[j]][3,]
  resplist_4[j,]=resplist[[j]][4,]
}
resplist_r = list()
resplist_r[[1]]=resplist_1
resplist_r[[2]]=resplist_2
resplist_r[[3]]=resplist_3
resplist_r[[4]]=resplist_4
return(resplist_r)
}