#######################################################################################################################
##  Adaptive Designs for Identifying Optimal Biological Dose for Molecularly Targeted Agents
##          by Yong Zang, Ying Yuan and Jack Lee
##
## This file contains six functions for implementation of the proposed logistic, isotonic and l-logistic designs to
## find the selected optimal biological dose (OBD) for targted agents. The OBD is defined as the lowest dose with the
## highest rate of efficacy while safe.
##
## logistic(), isotonic() and llogistic() are functions used to generate operating characteristics for the proposed designs
## df.logistic(), df.isotonic() and df.logistic() are functions used for dose-finding of the actual trial. 
##
## Based on numerical studies, we recommend isotonic() and llogistic() for use in practice because of their robust operating characteristics
## At the end of this file, an example is provided to illustrate how to use the proposed designs to conduct a clinical trial.
########################################################################################################################


##########################################################################################################
## Function to generate operating characteristics of the isotonic design 
## To use this function, library "Iso" should be installed in R. 
##
##  Arguments:
## ttox : toxicity rates for each dose level
## teff: efficacy rate for each dose level
## cohortsize: sample size for each cohort
## ncohort: total number of cohort
## ntrial: the number of simulated trial
## phi: upper bound of toxicity rate
## ct: threshold for posterior probability of toxicity; dose with toxicity probability larger than ct is excluded from the admissible set
#########################################################################################################

isotonic=function(ttox,teff,cohortsize=3,ncohort=10,ntrial=5000,phi=0.3,ct=0.8){
	library(Iso)
	
## find alpha in beta distribution
	f=function(alpha){
		a=pbeta(0.3,alpha,0.5-alpha)
		return(a-0.22)
	}
	alpha=uniroot(f,c(0.01,0.49))[[1]]
	
## find the admissible set
	findadm=function(n,ytox){
		a=alpha+ytox
		b=0.5-alpha+n-ytox
		re=1-pbeta(phi,a,b)
		re=pava(re)
		return(max(1,length(re[re<=ct])))
	}
	ndose=length(ttox)
	N=matrix(rep(0,ndose*ntrial),ncol=ndose)
	YTOX=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	YEFF=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	dselect=rep(0,ntrial)
	for (trial in 1: ntrial){
		ytox=rep(0,ndose)
		yeff=rep(0,ndose)
		n=rep(0,ndose)
		d=1 ##dose start from level 1
		adm=findadm(n,ytox)
		
## dose-finding procedure
		for (i in 1: ncohort){
			ytox[d]=ytox[d]+rbinom(1,cohortsize,ttox[d])
			yeff[d]=yeff[d]+rbinom(1,cohortsize,teff[d])
			n[d]=n[d]+cohortsize
			if (d==ndose){
				zfit=yeff/n
				unifit=ufit(zfit,x=c(1:ndose),type="b")[[2]]
				nd=tail(which(unifit==max(unifit)),1)
				if (nd==ndose){d=min(adm,d)}
				else {d=min(adm,d-1)}
			} else {
				try=length(n[n>0])
				zfit=yeff[1:try]/n[1:try]
				unifit=ufit(zfit,x=c(1:try),type="b")[[2]]
				nd=tail(which(unifit==max(unifit)),1)
				if ((nd==try)||(nd>d)){d=min(adm,d+1)}
				else if (nd<d) {d=min(adm,d-1)}
				else {d=min(adm,d)}
			}   
			adm=findadm(n,ytox)
		}
		z=yeff/(n+0.0001)
		unimode=ufit(z,x=c(1:ndose),type="b")[[2]]
		fid=which(unimode==max(unimode))
		dselect[trial]=min(fid[length(fid)],adm)
		N[trial,]=n
		YTOX[trial,]=ytox
		YEFF[trial,]=yeff
	}
	selpercent=rep(0, ndose)
	
## Summarize results 
	for(i in 1:ndose) { selpercent[i]=sum(dselect==i)/ntrial*100 }
	print("selection probablity") 
	cat(formatC(selpercent, digits=1, format="f"), sep="  ", "\n")
	print("average number of patients") 
	cat(formatC(c(apply(N,2,mean),sum(apply(N,2,mean))), digits=1, format="f"), sep ="  ", "\n")  
	print("average number of patients response to efficacy")
	cat(formatC(c(apply(YEFF,2,mean),sum(apply(YEFF,2,mean))), digits=1, format="f"), sep ="  ", "\n")
	print("average number of patients response to toxicity")
	cat(formatC(c(apply(YTOX,2,mean),sum(apply(YTOX,2,mean))), digits=1, format="f"), sep ="  ", "\n")
	
}

############################################################################
## Function for real trial dose-finding procedure of isotonic design
##
## n: number of patients treated at each dose level
## ytox: number of patients reporting toxicity at each dose level
## ytox: number of patients reporting efficacy at each dose level
############################################################################

df.isotonic=function(n,ytox,yeff,d,phi=0.3,ct=0.8){
	library(Iso)
## find alpha in beta distribution
	
	f=function(alpha){
		a=pbeta(0.3,alpha,0.5-alpha)
		return(a-0.22)
	}
	
	alpha=uniroot(f,c(0.01,0.49))[[1]]
	
## find the admissible set
	
	findadm=function(n,ytox){
		a=alpha+ytox
		b=0.5-alpha+n-ytox
		re=1-pbeta(phi,a,b)
		re=pava(re)
		return(max(1,length(re[re<=ct])))
	}
	ndose=length(ytox)
	adm=findadm(n,ytox)
	if (d==ndose){
		zfit=yeff/n
		unifit=ufit(zfit,x=c(1:ndose),type="b")[[2]]
		nd=tail(which(unifit==max(unifit)),1)
		if (nd==ndose){d=min(adm,d)}
		else {d=min(adm,d-1)}
	} else {
		try=length(n[n>0])
		zfit=yeff[1:try]/n[1:try]
		unifit=ufit(zfit,x=c(1:try),type="b")[[2]]
		nd=tail(which(unifit==max(unifit)),1)
		if ((nd==try)||(nd>d)){d=min(adm,d+1)}
		else if (nd<d) {d=min(adm,d-1)}
		else {d=min(adm,d)}
	}   
	
	return(list("next dose"=d))
	
}


##########################################################################################################
## Function to generate operating characteristics of the l-logistic design 
##
##  Arguments:
## ttox : toxicity rates for each dose level
## teff: efficacy rate for each dose level
## cohortsize: sample size for each cohort
## ncohort: total number of cohort
## ntrial: the number of simulated trial
## phi: upper bound of toxicity rate
## ct: threshold for posterior probability of toxicity; dose with toxicity probability larger than ct is excluded from the admissible set
## ce1, ce2: thresholds for the posterior probability of the slop
##           If the posterior probability of the slop is less than ce1, de-escalate the dose level
##           If the posterior probability of the slop is larger than ce2, escalate the dose level
##           Otherwise, retain at the current dose level
#########################################################################################################


llogistic=function(ttox,teff,cohortsize=3,ncohort=10,ntrial=5000,phi=0.3,ct=0.8,ce1=0.3,ce2=0.4){
	library(mcmc)
	
	library(Iso)
	
## find alpha in beta distribution
	
	f=function(alpha){
		a=pbeta(0.3,alpha,0.5-alpha)
		return(a-0.22)
	}
	
	alpha=uniroot(f,c(0.01,0.49))[[1]]
	
## find the admissible set
	
	findadm=function(n,ytox){
		a=alpha+ytox
		b=0.5-alpha+n-ytox
		re=1-pbeta(phi,a,b)
		re=pava(re)
		return(max(1,length(re[re<=ct])))
	}
	
## joint distribution for mcmc (efficacy)
	lupost <- function(beta, x, y) {
		eta <- cbind(1,x)%*%beta
		p <- 1/(1 + exp(-eta))
		logl <- sum(log(p[y == 1])) + sum(log(1-p[y == 0]))
		return(logl + dcauchy(beta[1], scale=10, log = TRUE)+dcauchy(beta[2], scale=2.5, log = TRUE))
	}
	
	collapse=function(a,b){
		re=NULL
		for (i in 1: length(a)){
			d=rep(c(1,0),c(a[i],b[i]-a[i]))
			re=c(re,d)
		}
		return(re)
	}
	thetaf=ce1
	thetab1=1-ce2
	thetab2=1-ce1
	ndose=length(ttox) ## number of dose levels
	x=c(1:ndose)
	xs=(x-mean(x))/(2*sd(x))
	N=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	YTOX=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	YEFF=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	dselect=rep(0,ntrial)
	step=ncohort-1
## start the dose finding
	for (trial in 1: ntrial){
		ytox=rep(0,ndose)
		yeff=rep(0,ndose)
		n=rep(0,ndose)
## put the first cohort on dose 1
		ytox[1]=rbinom(1,cohortsize,ttox[1])
		yeff[1]=rbinom(1,cohortsize,teff[1])
		n[1]=cohortsize
		adm=findadm(n,ytox) ## find the admissible set
		d=2 ## start for dose level 2
## dose-finding procedure
		for (i in 1: step){
			
			if (adm==1){
				d=1
				ytox[d]=ytox[d]+rbinom(1,cohortsize,ttox[d])
				yeff[d]=yeff[d]+rbinom(1,cohortsize,teff[d])
				n[d]=n[d]+cohortsize
				adm=findadm(n,ytox)
			} else {
				
				ytox[d]=ytox[d]+rbinom(1,cohortsize,ttox[d])
				yeff[d]=yeff[d]+rbinom(1,cohortsize,teff[d])
				n[d]=n[d]+cohortsize
				if (d==1){
					yfreg=collapse(yeff[d:(d+1)],n[d:(d+1)])
					xfreg=rep(xs[d:(d+1)],n[d:(d+1)])
					mcf=metrop(lupost,c(0,0),2000,x=xfreg,y=yfreg)
					mcf=mcf$batch[1001:2000,]
					if (mean(mcf[,2]>0)>thetaf){d=min((d+1),adm)}
					else (d=min(d,adm))
				} else if (d==ndose){
					ybreg=collapse(yeff[(d-1):d],n[(d-1):d])
					xbreg=rep(xs[(d-1):d],n[(d-1):d])
					mcb=metrop(lupost,c(0,0),2000,x=xbreg,y=ybreg)
					mcb=mcb$batch[1001:2000,]
					if (mean(mcb[,2]<=0)>thetab2){d=min((d-1),adm)}
					else (d=min(d,adm))
				} else if (n[d+1]==0){
					ybreg=collapse(yeff[(d-1):d],n[(d-1):d])
					xbreg=rep(xs[(d-1):d],n[(d-1):d])
					mcb=metrop(lupost,c(0,0),2000,x=xbreg,y=ybreg)
					mcb=mcb$batch[1001:2000,]
					if (mean(mcb[,2]<=0)>thetab2){d=min(d-1,adm)}
					else if (mean(mcb[,2]<=0)<=thetab1) (d=min(d+1,adm))
					else (d=min(d,adm))
				} else {
					yfreg=collapse(yeff[d:(d+1)],n[d:(d+1)])
					xfreg=rep(xs[d:(d+1)],n[d:(d+1)])
					mcf=metrop(lupost,c(0,0),2000,x=xfreg,y=yfreg)
					mcf=mcf$batch[1001:2000,]
					ybreg=collapse(yeff[(d-1):d],n[(d-1):d])
					xbreg=rep(xs[(d-1):d],n[(d-1):d])
					mcb=metrop(lupost,c(0,0),2000,x=xbreg,y=ybreg)
					mcb=mcb$batch[1001:2000,]
					if ((mean(mcb[,2]<=0)<=thetab1)&&(mean(mcf[,2]>0)>thetaf)){d=min(d+1,adm)}
					else if (mean(mcb[,2]<=0)>thetab2){d=min(d-1,adm)}
					else (d=min(d,adm))
					
				}
				adm=findadm(n,ytox)
			}
		}
		zeff=yeff/(n+0.0001)
		unimode=ufit(zeff,x=c(1:ndose),type="b")[[2]]
		fid=which(unimode==max(unimode))
		dselect[trial]=min(fid[length(fid)],adm)
		N[trial,]=n
		YTOX[trial,]=ytox
		YEFF[trial,]=yeff
	}
## record the results
	selpercent=rep(0, ndose)
	for(i in 1:ndose) { selpercent[i]=sum(dselect==i)/ntrial*100 }
	print("selection probablity") 
	cat(formatC(selpercent, digits=1, format="f"), sep="  ", "\n")
	print("average number of patients") 
	cat(formatC(c(apply(N,2,mean),sum(apply(N,2,mean))), digits=1, format="f"), sep ="  ", "\n")
	print("average number of patients response to efficacy")
	cat(formatC(c(apply(YEFF,2,mean),sum(apply(YEFF,2,mean))), digits=1, format="f"), sep ="  ", "\n")
	print("average number of patients response to toxicity")
	cat(formatC(c(apply(YTOX,2,mean),sum(apply(YTOX,2,mean))), digits=1, format="f"), sep ="  ", "\n")
	
}

############################################################################
## Function for real trial dose-finding procedure of l-logistic design
#
## n: number of patients treated at each dose level
## ytox: number of patients reporting toxicity at each dose level
## ytox: number of patients reporting efficacy at each dose level
############################################################################


df.llogistic=function(n,ytox,yeff,d,phi=0.3,ct=0.8,ce1=0.3,ce2=0.4){
	library(mcmc)
	
	library(Iso)
	
## find alpha in beta distribution
	
	f=function(alpha){
		a=pbeta(0.3,alpha,0.5-alpha)
		return(a-0.22)
	}
	
	alpha=uniroot(f,c(0.01,0.49))[[1]]
	
## find the admissible set
	
	findadm=function(n,ytox){
		a=alpha+ytox
		b=0.5-alpha+n-ytox
		re=1-pbeta(phi,a,b)
		re=pava(re)
		return(max(1,length(re[re<=ct])))
	}
	
## joint distribution for mcmc (efficacy)
	lupost <- function(beta, x, y) {
		eta <- cbind(1,x)%*%beta
		p <- 1/(1 + exp(-eta))
		logl <- sum(log(p[y == 1])) + sum(log(1-p[y == 0]))
		return(logl + dcauchy(beta[1], scale=10, log = TRUE)+dcauchy(beta[2], scale=2.5, log = TRUE))
	}
	
	collapse=function(a,b){
		re=NULL
		for (i in 1: length(a)){
			d=rep(c(1,0),c(a[i],b[i]-a[i]))
			re=c(re,d)
		}
		return(re)
	}
	thetaf=ce1
	thetab1=1-ce2
	thetab2=1-ce1
	ndose=length(ytox) ## number of dose levels
	x=c(1:ndose)
	xs=(x-mean(x))/(2*sd(x))
	adm=findadm(n,ytox)
	if (adm==1){
		d=1
	} else {  
		if (d==1){
			yfreg=collapse(yeff[d:(d+1)],n[d:(d+1)])
			xfreg=rep(xs[d:(d+1)],n[d:(d+1)])
			mcf=metrop(lupost,c(0,0),2000,x=xfreg,y=yfreg)
			mcf=mcf$batch[1001:2000,]
			if (mean(mcf[,2]>0)>thetaf){d=min((d+1),adm)}
			else (d=min(d,adm))
		} else if (d==ndose){
			ybreg=collapse(yeff[(d-1):d],n[(d-1):d])
			xbreg=rep(xs[(d-1):d],n[(d-1):d])
			mcb=metrop(lupost,c(0,0),2000,x=xbreg,y=ybreg)
			mcb=mcb$batch[1001:2000,]
			if (mean(mcb[,2]<=0)>thetab2){d=min((d-1),adm)}
			else (d=min(d,adm))
		} else if (n[d+1]==0){
			ybreg=collapse(yeff[(d-1):d],n[(d-1):d])
			xbreg=rep(xs[(d-1):d],n[(d-1):d])
			mcb=metrop(lupost,c(0,0),2000,x=xbreg,y=ybreg)
			mcb=mcb$batch[1001:2000,]
			if (mean(mcb[,2]<=0)>thetab2){d=min(d-1,adm)}
			else if (mean(mcb[,2]<=0)<=thetab1) (d=min(d+1,adm))
			else (d=min(d,adm))
		} else {
			yfreg=collapse(yeff[d:(d+1)],n[d:(d+1)])
			xfreg=rep(xs[d:(d+1)],n[d:(d+1)])
			mcf=metrop(lupost,c(0,0),2000,x=xfreg,y=yfreg)
			mcf=mcf$batch[1001:2000,]
			ybreg=collapse(yeff[(d-1):d],n[(d-1):d])
			xbreg=rep(xs[(d-1):d],n[(d-1):d])
			mcb=metrop(lupost,c(0,0),2000,x=xbreg,y=ybreg)
			mcb=mcb$batch[1001:2000,]
			if ((mean(mcb[,2]<=0)<=thetab1)&&(mean(mcf[,2]>0)>thetaf)){d=min(d+1,adm)}
			else if (mean(mcb[,2]<=0)>thetab2){d=min(d-1,adm)}
			else (d=min(d,adm))
			
		}
		
	}
	return(list("next dose"=d))
}




###########################################################################################################################################
## Function to generate operating characteristics of the logistic design 
## ttox : toxicity rates for each dose level
## teff: efficacy rate for each dose level
## cohortsize: sample size for each cohort
## ncohort: total number of cohort
## ntrial: the number of simulated trial
## phi: upper bound of toxicity rate
## ct: threshold for posterior probability of toxicity; dose with toxicity probability larger than ct is excluded from the admissible set
##############################################################################################################################################

logistic=function(ttox,teff,cohortsize=3,ncohort=10,ntrial=5000,phi=0.3,ct=0.8){
	library(arm)  
	library(Iso)
	
## find alpha in beta distribution
	
	f=function(alpha){
		a=pbeta(0.3,alpha,0.5-alpha)
		return(a-0.22)
	}
	
	alpha=uniroot(f,c(0.01,0.49))[[1]]
	
## find the admissible set
	
	findadm=function(n,ytox){
		a=alpha+ytox
		b=0.5-alpha+n-ytox
		re=1-pbeta(phi,a,b)
		re=pava(re)
		return(max(1,length(re[re<=ct])))
	}
	
## convert binomial to bernoulli
	
	collapse=function(a,b){
		re=NULL
		for (i in 1: length(a)){
			d=rep(c(1,0),c(a[i],b[i]-a[i]))
			re=c(re,d)
		}
		return(re)
	}
	
	ndose=length(ttox) 
	x=c(1:ndose)
	x1s=(x-mean(x))/(2*sd(x))
	x2s=x1s^2
	N=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	YTOX=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	YEFF=matrix(rep(0,ndose*ntrial),ncol=ndose) 
	dselect=rep(0,ntrial) 
	for (trial in 1:ntrial){
		ytox=rep(0,ndose)
		yeff=rep(0,ndose)
		d=1
		n=rep(0,ndose)
		adm=findadm(n,ytox)
## dose-finding procedure ############
		for (i in 1: ncohort){
			ytox[d]=ytox[d]+rbinom(1,cohortsize,ttox[d])
			yeff[d]=yeff[d]+rbinom(1,cohortsize,teff[d])
			n[d]=n[d]+cohortsize
			yreg=collapse(yeff,n)
			x1reg=rep(x1s,n)
			x2reg=rep(x2s,n)
			mc=bayesglm(yreg~x1reg+x2reg,family=binomial(link="logit")) 
			est=1/(1+exp(-mc[[1]][1]-mc[[1]][2]*x1s-mc[[1]][3]*x2s))
			nd=which(est==max(est))
			nd=nd[length(nd)]
			if(nd>d){d=min((d+1),ndose,adm)}
			else if(nd<d) {d=min(adm,max(1,(d-1)))}
			else (d=min(adm,d))
			adm=findadm(n,ytox)
		}
		dselect[trial]=min(nd,adm)
		N[trial,]=n
		YTOX[trial,]=ytox
		YEFF[trial,]=yeff
	}
## Summarize results
	selpercent=rep(0, ndose)
	for(i in 1:ndose) { selpercent[i]=sum(dselect==i)/ntrial*100 }
	print("selection probablity") 
	cat(formatC(selpercent, digits=1, format="f"), sep="  ", "\n")
	print("average number of patients") 
	cat(formatC(c(apply(N,2,mean),sum(apply(N,2,mean))), digits=1, format="f"), sep ="  ", "\n")  
	print("average number of patients response to efficacy")
	cat(formatC(c(apply(YEFF,2,mean),sum(apply(YEFF,2,mean))), digits=1, format="f"), sep ="  ", "\n")
	print("average number of patients response to toxicity")
	cat(formatC(c(apply(YTOX,2,mean),sum(apply(YTOX,2,mean))), digits=1, format="f"), sep ="  ", "\n")
	
}

############################################################################
## Function for real trial dose-finding procedure of logistic design
## n: number of patients treated at each dose level
## ytox: number of patients reporting toxicity at each dose level
## ytox: number of patients reporting efficacy at each dose level
###############################################################################
df.logistic=function(n,ytox,yeff,d,phi=0.3,ct=0.8){
	library(arm)  
	library(Iso)
	
## find alpha in beta distribution
	
	f=function(alpha){
		a=pbeta(0.3,alpha,0.5-alpha)
		return(a-0.22)
	}
	
	alpha=uniroot(f,c(0.01,0.49))[[1]]
	
## find the admissible set
	
	findadm=function(n,ytox){
		a=alpha+ytox
		b=0.5-alpha+n-ytox
		re=1-pbeta(phi,a,b)
		re=pava(re)
		return(max(1,length(re[re<=ct])))
	}
	
## convert binomial to bernoulli
	
	collapse=function(a,b){
		re=NULL
		for (i in 1: length(a)){
			d=rep(c(1,0),c(a[i],b[i]-a[i]))
			re=c(re,d)
		}
		return(re)
	}
	ndose=length(ytox) ## number of total dose levels 
	x=c(1:ndose)
	x1s=(x-mean(x))/(2*sd(x))
	x2s=x1s^2
	adm=findadm(n,ytox)
	yreg=collapse(yeff,n)
	x1reg=rep(x1s,n)
	x2reg=rep(x2s,n)
	mc=bayesglm(yreg~x1reg+x2reg,family=binomial(link="logit")) 
	est=1/(1+exp(-mc[[1]][1]-mc[[1]][2]*x1s-mc[[1]][3]*x2s))
	nd=which(est==max(est))
	nd=nd[length(nd)]
	if(nd>d){d=min((d+1),ndose,adm)}
	else if(nd<d) {d=min(adm,max(1,(d-1)))}
	else (d=min(adm,d))
	
	return(list("next dose"=d))
	
}




########################################## an example #######################################

####### generate operating characteristics of the designs ##########
# assume that the target toxicity rate is 30%, the sample size is 30 in cohorts of size 3
eff=c(0.05,0.25,0.6,0.8,0.5) ## true efficacy rate
tox=c(0.1,0.2,0.45,0.55,0.65) ## true toxicity rate 
isotonic(ttox=tox,teff=eff,cohortsize=3,ncohort=10,ntrial=50,phi=0.3,ct=0.8) ##  generate operating characteristics of isotoinc design
llogistic(ttox=tox,teff=eff,cohortsize=3,ncohort=10,ntrial=50,phi=0.3,ct=0.8,ce1=0.3,ce2=0.4) ##  generate operating characteristics of l-logistic design



#############  trial conduct  #############
# assume that the first cohort is treated at dose level 1 without toxicity and one efficacy response, yielding the data
n=c(3,0,0,0,0)
ytox=c(0,0,0,0,0)
yeff=c(1,0,0,0,0)
d=1

## call dose-finding function to obtain the recommmended dose for the next cohort
df.isotonic(n,ytox,yeff,d,phi=0.3,ct=0.8)   ## using the isotonic method
##$`next dose`
##[1] 2
df.llogistic(n,ytox,yeff,d,phi=0.3,ct=0.8,ce1=0.3,ce2=0.4)   ## or use the l-logistic method
##$`next dose`
##[1] 2

## Hence, we escalate the dose level and treat the second cohort at dose level 2
## assume that in cohort 2, no patient experience toxicity and 2 patients have efficacy response at dose level 2, yielding the updated data
n=c(3,3,0,0,0)
ytox=c(0,0,0,0,0)
yeff=c(1,2,0,0,0)
d=2
## call dose-finding function to obtain the recommmended dose for the next cohort
df.isotonic(n,ytox,yeff,d,phi=0.3,ct=0.8) ## using the isotonic method
##$`next dose`
##[1] 3
df.llogistic(n,ytox,yeff,d,phi=0.3,ct=0.8,ce1=0.3,ce2=0.4)  ## or use the l-logistic method
##$`next dose`
##[1] 3

## Hence, we escalate the dose level and treat the thrid cohort at dose level 3
.
.
## ....... repeat the above procedure for incoming cohorts................
.
.
## Suppose after the last cohort of patients have been treated, the updated data are
n=c(3,6,12,6,3)
ytox=c(0,0,1,1,2)
yeff=c(1,2,6,1,0)
d=3

## call dose-finding function to obtain the recommmended optimal biological dose (OBD).
df.isotonic(n,ytox,yeff,d,phi=0.3,ct=0.8) 
##$`next dose`
##[1] 3
df.llogistic(n,ytox,yeff,d,phi=0.3,ct=0.8,ce1=0.3,ce2=0.4)  
##$`next dose`
##[1] 3

## Both the isotonic and l-logistic methods select dose 3 as the OBD.