##### The R code to implement the Bayesian hybrid design by Yuan and Yin (Statistics in Medicine, 2010) ##### ################################### parameters for the function ################################# # target: the target toxicity probability # delta: the tolerable margin for the MTD # ptox: true toxicity probabilities # cohortsize: the size of cohort # ncohort: the number of cohorts for the trial # ntrial: the number of simulated trials # p: the skeleton for the CRM model, that is, Pr(toxicity at dose j) = p_j^(exp(alpha)) ################################################################################################## hybridf <- function(target, delta, ptox, cohortsize, ncohort, ntrial=1000, p) { ##################### subroutines #################### # marginal density marginalf <- function(alpha, p, y, n, d) { lik=1; for(j in 1:length(p)) { pj = p[j]^(exp(alpha)); lik = lik*pj^y[j]*(1-pj)^(n[j]-y[j]); } return(lik*exp(alpha)*p[d]^(exp(alpha))); } # posterior = likelihood x prior posterior <- function(alpha, p, y, n) { sigma2 = 2; lik=1; for(j in 1:length(p)) { pj = p[j]^(exp(alpha)); lik = lik*pj^y[j]*(1-pj)^(n[j]-y[j]); } return(lik*exp(-0.5*alpha*alpha/sigma2)); } # the posterior mean of ptox posttoxf <- function(alpha, p, y, n, j) { p[j]^(exp(alpha))*posterior(alpha, p, y, n); } # isotonic transformation using the pool adjacent violator algorithm (PAVA) 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 } ##################### end of subroutines #################### ndose=length(ptox); Y=matrix(rep(0, ndose*ntrial), ncol=ndose); N=matrix(rep(0, ndose*ntrial), ncol=ndose); ptox.hat = numeric(ndose); # estimate of toxicity prob dselect = rep(0, ntrial); TOXSTOP = 0.9; ## toxicity stop boundary when using CRM rule lambda = 0.61; ################## simulate trials ################### for(trial in 1:ntrial) { y<-rep(0,ndose); n<-rep(0,ndose); earlystop=0; d=1; # dose level starting from 1 maxdose=1; # highest for(i in 1:ncohort) { ### generate random toxicity response y[d] = y[d] + sum(runif(cohortsize)TOXSTOP) { earlystop=1; break;} ### calculate bayes factors mlik1 = pbeta(target-delta, y[d]+1, n[d]-y[d]+1)/(target-delta); # H1: toxicitytarget+delta pH1 = mlik1/(mlik1+mlik2+mlik3); pH2 = mlik2/(mlik1+mlik2+mlik3); pH3 = mlik3/(mlik1+mlik2+mlik3); ### dose escalation/de-escalation if(pH1>lambda && d!=ndose) { d = d+1; } else if(pH3>lambda && d!=1) { d = d-1; } else if(pH2>lambda) { d = d; } else ## switch to CRM model-based approach for dose assignment { ## calculate model-based BF b1 = log(log(target-delta)/log(p[d])) b2 = log(log(target+delta)/log(p[d])) cmlik1 = integrate(marginalf,lower=b1,upper=10,p,y,n,d)$value/(target-delta); cmlik2 = integrate(marginalf,lower=b2,upper=b1,p,y,n,d)$value/(2*delta) cmlik3 = integrate(marginalf,lower=-Inf,upper=b2,p,y,n,d)$value/(1-target-delta); cpH1 = cmlik1/(cmlik1+cmlik2+cmlik3); cpH2 = cmlik2/(cmlik1+cmlik2+cmlik3); cpH3 = cmlik3/(cmlik1+cmlik2+cmlik3); # dose escalation/de-escalation if(cpH1>lambda && d!=ndose) { d=d+1; } else if(cpH3>lambda && d!=1) { d=d-1; } else if(cpH2>lambda) { d=d; } } maxdose=max(maxdose, d); # recording the highest tried dose } Y[trial,]=y; N[trial,]=n; if(earlystop==1) { dselect[trial]=99; } else ## selection the MTD based on isotonic transformation { pp<-rep(0, maxdose) pp.var<-rep(0, maxdose) for(i in 1:maxdose){ # using the prior Beta(0.005, 0.005) to obtain posterior estimates pp[i] <- (y[i]+.005)/(n[i]+.01); pp.var[i] <- (y[i]+0.005)*(n[i]-y[i]+0.005)/((n[i]+0.01)^2*(n[i]+0.01+1)) } pp<-pava(pp, wt=1/pp.var) ## perform the isotonic transformation using PAVA for(i in 2:maxdose){ pp[i] <- pp[i] + i*1E-10} ## break ties by adding an increasingly small number seldose<-sort(abs(pp-target), index.return=T)$ix[1] ## select dose closest to the target as the MTD dselect[trial] <- seldose; ## recoding dose selection } } ################## output results ################################ selpercent=rep(0, ndose); for(i in 1:ndose) { selpercent[i]=sum(dselect==i)/ntrial*100; } # selection probablity cat(formatC(c(selpercent, sum(dselect==99)/ntrial*100), digits=1, format="f"), sep=" ", "\n"); # print the average number of patients cat(formatC(c(apply(N,2,mean), sum(N[,ptox>target])/ntrial, sum(Y)/ntrial), digits=1, format="f"), sep =" ", "\n"); } #### An example: scenario 1 in the paper #### ncohort=8; cohortsize=3; target=0.30; delta=target/10; p=c(0.14, 0.20, 0.25, 0.30, 0.35, 0.40); ptox <- c(0.10, 0.12, 0.30, 0.50, 0.60, 0.65) hybridf(target, delta, ptox, cohortsize, ncohort, ntrial=1000, p)