######################################################################### # # This R Program simulates the method of Wages and Tait (2015)* # # *Wages, N. and Tait, C. (2015) Seamless Phase I/II Adaptive Design # for Oncology Trials of Molecularly Targeted Therapies. J Biopharm Stats. # ######################################################################### ###install required R packages library(nnet) library(dfcrm) library(binom) ###Load the function 'bpocrm' bpocrm<-function(p0,q0,p.skel,q.skel,tul,ell,cohortsize,ncohort,start.comb){ # if a single ordering is inputed as a vector, convert it to a matrix if(is.vector(q.skel)) q.skel=t(as.matrix(q.skel)); nord.eff = nrow(q.skel); mprior.eff = rep(1/nord.eff, nord.eff); # prior for each efficacy ordering avar=1.34 bcrmh<-function(a,p,y,n){ s2=avar lik=exp(-0.5*a*a/s2) for(j in 1:length(p)){ pj=p[j]**exp(a) lik=lik*pj^y[j]*(1-pj)^(n[j]-y[j]); } return(lik); } bcrmht<-function(a,p,y,n){ s2=avar lik=a*exp(-0.5*a*a/s2) for(j in 1:length(p)){ pj=p[j]**exp(a) lik=lik*pj^y[j]*(1-pj)^(n[j]-y[j]); } return(lik); } bcrmht2<-function(a,p,y,n){ s2=avar lik=a^2*exp(-0.5*a*a/s2) for(j in 1:length(p)){ pj=p[j]**exp(a) lik=lik*pj^y[j]*(1-pj)^(n[j]-y[j]); } return(lik); } ### run a trial ncomb = ncol(q.skel); #number of combos y=rep(0,ncomb); #number of toxicity/responses at each dose level z=rep(0,ncomb); #number of efficacy at each dose level n=rep(0,ncomb); #number of treated patients at each dose level comb.curr = start.comb; # current dose level ptox.hat = numeric(ncomb); # estimate of toxicity prob comb.select=rep(0,ncomb); # a vector of indicators for dose selection stop=stops=stopf=0; #indicate if trial stops early i=1 while(i <= ncohort) { # generate data for a new cohort of patients y[comb.curr] = y[comb.curr] + rbinom(1,cohortsize,p0[comb.curr]); z[comb.curr] = z[comb.curr] + rbinom(1,cohortsize,q0[comb.curr]); n[comb.curr] = n[comb.curr] + cohortsize; ####Toxicity marginal= integrate(bcrmh,lower=-Inf,upper=Inf, p=p.skel, y=y,n=n,abs.tol = 0)\$value; est=integrate(bcrmht,lower=-10,upper=10, p.skel, y, n,abs.tol = 0)\$value/marginal ptox.hat=p.skel**exp(est) ####Efficacy marginal.eff = est.eff=rep(0, nord.eff); for(k in 1:nord.eff) { marginal.eff[k] = integrate(bcrmh,lower=-Inf,upper=Inf, p=q.skel[k,], y=z,n=n,abs.tol = 0)\$value; est.eff[k]=integrate(bcrmht,lower=-10,upper=10, q.skel[k,], z, n,abs.tol = 0)\$value/marginal.eff[k] } postprob.eff = (marginal.eff*mprior.eff)/sum(marginal.eff*mprior.eff); # efficacy model selection, identify the model with the highest posterior prob if(nord.eff>1){ meff.sel = which.is.max(postprob.eff); } else{ meff.sel = 1; } aset=which(ptox.hat<=tul) if(length(aset)==0){aset=which.min(ptox.hat)} peff.hat=matrix(,nrow=nrow(q.skel),ncol=d) peff.hat=q.skel**exp(est.eff) peff.hat.aset=rep(0,ncomb) peff.hat.aset[aset]=peff.hat[meff.sel,aset] ri0=peff.hat.aset/sum(peff.hat.aset) if(length(aset)==1){ comb.best=aset } else { ifelse(sum(n)<=n.ar,comb.best<-sample(1:ncomb,1,prob=ri0),comb.best<-which.max(peff.hat.aset)) } ##########stopping rules safety=binom.confint(y[1],n[1],conf.level=0.95,methods="exact")\$lower if(safety>tul){ stops=1 break } if(sum(n) > n.ar){ futility=binom.confint(z[comb.best],n[comb.best],conf.level=0.95,methods="exact")\$upper if(futility