####################################################################################
#									
# 	This R Program simulates the design of Horton, Wages, Gentzler (2021+)* 				
#									
#	*Horton BJ, Wages NA, Gentzler RD. Bayesian early-phase design
#	 for identifying cohort-specific optimal dose combinations based on 
#	 multiple endpoints: application to a phase I trial in Non-Small Cell
#	 Lung Cancer. Int J Environ Res Public Health 2021+.						
#									
###################################################################################

###install required R packages
#install.packages("nnet")
#install.packages("dfcrm")
library(nnet)
library(dfcrm)


###Load the function 'HWGdesign' 
HWGdesign<-function(p0,q0,p.skel,tul,ell,cohortsize,ncohort,n.stop,start.comb,cs,cf,avar){

n.ar<-ncohort*cohortsize/3 

# if a single ordering is inputed as a vector, convert it to a matrix
	if(is.vector(p.skel)) p.skel=t(as.matrix(p.skel));
	
	nord.tox = nrow(p.skel);
	mprior.tox = rep(1/nord.tox, nord.tox);  # prior for each toxicity ordering

post.tox<-function(a,p,y,n){
	s2=avar
	lik=1
	for(j in 1:length(p)){
		pj=p[j]**exp(a)
		lik=lik*pj^y[j]*(1-pj)^(n[j]-y[j]);
		}
	return(lik*exp(-0.5*a*a/s2));
    }

#the posterior mean of ptox
    posttoxf <- function(a, p, y, n, j) { p[j]^(exp(a))*post.tox(a, p, y, n); }

### run a trial 	
    ncomb = ncol(p.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
    peff.hat = numeric(ncomb); # estimate of efficacy prob
    p.safety = numeric(ncomb)
    p.futility = numeric(ncomb)
    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)
    {

		if(n[comb.curr]==n.stop){
			stop<-0
			break
		}
	# generate data for a new cohort of patients
		lasttox=rbinom(1,cohortsize,p0[comb.curr]);
		y[comb.curr] = y[comb.curr] + lasttox;
		z[comb.curr] = z[comb.curr] + rbinom(1,cohortsize,q0[comb.curr]);
		n[comb.curr] = n[comb.curr] + cohortsize;		

		

		marginal.tox = rep(0, nord.tox);
		for(k in 1:nord.tox)
		{
			marginal.tox[k] = integrate(post.tox,lower=-Inf,upper=Inf, p=p.skel[k,], y=y,n=n)$value;
		}
		
		postprob.tox = (marginal.tox*mprior.tox)/sum(marginal.tox*mprior.tox);

		
		
		# toxicity model selection, identify the model with the highest posterior prob
			if(nord.tox>1){ 
				mtox.sel = which.max(postprob.tox); 
			} else{
				mtox.sel = 1;
			}

		

		# calculate posterior mean of toxicity probability at each combo
			for(j in 1:ncomb){
				ptox.hat[j] = integrate(posttoxf,lower=-Inf,upper=Inf, p.skel[mtox.sel,], y, n,j)$value/marginal.tox[mtox.sel]; 
				p.safety[j] = integrate(post.tox,lower=log(log(tul)/log(p.skel[mtox.sel,j]+0.01)),upper=Inf, p.skel[mtox.sel,], y, n)$value/marginal.tox[mtox.sel];
			}
		

		peff.hat<-(z+0.5)/(n+1)
		
		if(p.safety[which.min(ptox.hat)]<cs){
				stops<-1
				break
				}

		mtd=which.is.max(-abs(ptox.hat-tul))
		utility=peff.hat*as.numeric(ptox.hat<=ptox.hat[mtd])
		cand=which(utility==max(utility))
		med=cand[ptox.hat[cand]==min(ptox.hat[cand])]
		ri=utility/sum(utility)

		comb.curr<-ifelse(sum(n)<n.ar,sample(1:ncomb,1,prob=ri),med)	

		if(sum(n)>n.ar){
			if(1 - pbeta(ell, z[comb.curr] + 0.5, n[comb.curr] - z[comb.curr] + 0.5) < cf){
				stopf<-1
				break
				}
			}		
	
		i=i+1
 		}

	if(stop==0 & stops==0 & stopf==0){
		
		comb.select[comb.curr]=comb.select[comb.curr]+1;
		}
	return(list(comb.select=comb.select,tox.data=y,eff.data=z,pt.allocation=n,stopf=stopf,stops=stops))
}
##########'HWGdesign' end here

###Load the function 'HWGdesign.sim' 
HWGdesign.sim<-function(pa0,qa0,ps0,qs0,p.skel,tul,ella,ells,cohortsize,ncohort,n.stop,ntrial,start.comb,cs,cf,pera,avar){
	ncomb=length(pa0)
		
	nptsa=round(ncohort*pera,0)
	nptss=ncohort-nptsa

	comb.selects<-ys<-zs<-ns<-matrix(nrow=ntrial,ncol=ncomb)
	nstopfs=nstopss=0

	comb.selecta<-ya<-za<-na<-matrix(nrow=ntrial,ncol=ncomb)
	nstopfa=nstopsa=0

	
	for(i in 1:ntrial){
		resulta<-HWGdesign(pa0,qa0,p.skel,tul,ella,cohortsize,nptsa,n.stop,start.comb,cs,cf,avar)
		comb.selecta[i,]=resulta$comb.select
		ya[i,]=resulta$tox.data
		za[i,]=resulta$eff.data
		na[i,]=resulta$pt.allocation
		nstopfa=nstopfa+resulta$stopf
		nstopsa=nstopsa+resulta$stops
		
		results<-HWGdesign(ps0,qs0,p.skel,tul,ells,cohortsize,nptss,n.stop,start.comb,cs,cf,avar)
		comb.selects[i,]=results$comb.select
		ys[i,]=results$tox.data
		zs[i,]=results$eff.data
		ns[i,]=results$pt.allocation
		nstopfs=nstopfs+results$stopf
		nstopss=nstopss+results$stops


	}

	cat("True toxicity probability in Adeno:\n")
	cat(round(pa0,3),  sep="\t",  "\n")
	cat("True activity probability in Adeno:\n")
	cat(round(qa0,3), sep="\t",  "\n")
	cat("Selection percentage in Adeno:\n")
	cat(formatC(colMeans(comb.selecta)*100, digits=1, format="f"), sep="\t",  "\n")
	cat("Average number of toxicities in Adeno:\n")
	cat(formatC(colMeans(ya), digits=1, format="f"), sep="\t",   "\n")
	cat("Average number of responses in Adeno:\n")
	cat(formatC(colMeans(za), digits=1, format="f"), sep="\t",   "\n")
	cat("Average number of patients treated in Adeno:\n")
	cat(formatC(colMeans(na), digits=1, format="f"), sep="\t",   "\n")
	cat("Percentage of trials stopped for safety in Adeno:\n")
	cat(nstopsa/ntrial*100, "\n")
	cat("Percentage of trials stopped for futility in Adeno:\n")
	cat(nstopfa/ntrial*100, "\n\n\n")
	cat("True toxicity probability in SCC:\n")
	cat(round(ps0,3),  sep="\t",  "\n")
	cat("True activity probability in SCC:\n")
	cat(round(qs0,3), sep="\t",  "\n")
	cat("Selection percentage in SCC:\n")
	cat(formatC(colMeans(comb.selects)*100, digits=1, format="f"), sep="\t",  "\n")
	cat("Average number of toxicities in SCC:\n")
	cat(formatC(colMeans(ys), digits=1, format="f"), sep="\t",   "\n")
	cat("Average number of responses in SCC:\n")
	cat(formatC(colMeans(zs), digits=1, format="f"), sep="\t",   "\n")
	cat("Average number of patients treated in SCC:\n")
	cat(formatC(colMeans(ns), digits=1, format="f"), sep="\t",   "\n")
	cat("Percentage of trials stopped for safety in SCC:\n")
	cat(nstopss/ntrial*100, "\n")
	cat("Percentage of trials stopped for futility in SCC:\n")
	cat(nstopfs/ntrial*100, "\n")

}
##########'HWGdesign.sim' end here




###############Trial simualtions
###number of study combinations
d<-6  

#####number of possible orderings
s<-4 ##for toxicity

###Specify the possible DLT orderings of the drug combinations
orders<-matrix(nrow=s,ncol=d)
orders[1,]<-c(1,4,2,5,3,6)	#down diagonals
orders[2,]<-c(1,2,4,3,5,6)	#up diagonals
orders[3,]<-c(1,4,2,3,5,6)	#alternate down-up diagonals
orders[4,]<-c(1,2,4,5,3,6)	#alternate down-up diagonals

###Specify skeleton for DLT probabilities
skeleton<-round(getprior(0.04,0.30,6,6),2)

###The following will adjust the location of each skeleton value to correspond 
###to the 's' possible orderings specified above. 
###'alpha' is the matrix of skeleton values
p.skel<-matrix(0,nrow=s,ncol=d)
for(j in 1:s){ 
	p.skel[j,]<-skeleton[order(orders[j,])]
}

###True probability scenarios
p1<-c(0.20,0.25,0.30,0.22,0.27,0.32)
p3<-c(0.01,0.05,0.07,0.08,0.10,0.15)
p2<-c(0.08,0.20,0.40,0.10,0.22,0.42)
p4<-c(0.01,0.05,0.15,0.03,0.08,0.20)

q1<-c(0.45,0.57,0.68,0.55,0.67,0.78)
q2<-c(0.55,0.83,0.68,0.60,0.85,0.70)
q3<-c(0.65,0.83,0.68,0.70,0.85,0.70)
q4<-c(0.25,0.40,0.35,0.35,0.50,0.45)


tul<-0.30     	   ##toxicity upper limit 
ella<-ells<-0.28     ##efficacy lower limit
cohortsize=1         ##cohort size for each inclusion
ncohort=60           ##number of cohorts
start.comb=2         ##starting combination
n.stop=12	         ##number needed on one combo to stop
ntrial=200           ##number of simulated trials 
cs=cf=0.10	         ##thresholds used to define acceptable combinations
pera=0.65		   ##anticipated participant accrual percentage from Cohort A 
avar <- 0.23	   ##variance of the dose-toxicity model parameter

##simulate a single trial
#HWGdesign(pa0,qa0,p.skel,tul,ella,cohortsize,nptsa,n.stop,start.comb,cs,cf,avar)

pa0<-p1
ps0<-p1
qa0<-q1
qs0<-q1
HWGdesign.sim(pa0,qa0,ps0,qs0,p.skel,tul,ella,ells,cohortsize,ncohort,n.stop,ntrial,start.comb,cs,cf,pera,avar)