#########################################################################
#									
# 	This R Program is for the simulation of the CRM 		
#	for partial orders in Phase I clinical trials with multiple treatment schedules.
#
#	The function 'crm' calculates MLE of parameter a and recommends a dose-schedule
#       combination for next patient based on power model.
#
#	The function 'twostgcrm' returns output from a single simulated trial. 	
#
#	The function 'pocrm.sim' returns the results of many 					
#	simulated trials.
#
#	Begin by loading these three functions.						
#							
#########################################################################

###LOAD FUNCTION 'crm'

crm<-function(obs,alpha,theta){

##'obs' is a data frame that conatins patient, dose and 
##toxicity information for those already in the study
##'alpha' is the matrix of alpha values for the working 
##model
##'theta' is the target MTC toxicity rate
##crm function yields the recommended dose for the next 
##patient

sim<<-table(obs$level,obs$tox)
ifelse(dim(sim)[2]==1,sim1<<-data.frame(as.numeric(row.names(sim)),sim[,1],1-sim[,1]),sim1<<-data.frame(as.numeric(row.names(sim)),sim[,1],sim[,2]))
names(sim1)<<-c('level','nontox','tox')

apred<<-rep(0,s)
	lik<<-rep(0,s)
	for(k in 1:s){			
		ll<-function(a){
			la<-0 #value of the log-likelihood
			for(i in sim1$level){
			index<-match(i,sim1$level)
			la<-la+sim1$tox[index]*a*log(alpha[k,][i])+sim1$nontox[index]*log((1-alpha[k,][i]**a))
			}
			la
			}
		apred[k]<<-optimize(f=ll,interval=c(0,100),maximum=T)$maximum
		lik[k]<<-ll(apred[k])
		}
pord<<-exp(lik)/sum(exp(lik))
library("nnet")
ord<<-which.is.max(pord)
ahat<<-apred[ord]
rpred<<-alpha[ord,]**ahat
next.lev<<-which.is.max(-(abs(rpred-theta)))
next.lev
    }
###'crm' ENDS HERE


###LOAD FUNCTION 'twostgcrm'

twostgcrm<-function(n,alpha,r,theta){
#returns the dose and tox info for each patient
#plus one more patient
#alpha is the vector of alpha values for the working 
#model
#r is the true toxicity rate for each dose
#theta is again the target MTD toxicity rate
#twostgcrm will generate the dose and toxicity 
#information for each patient, plus one more

	n1<<-n+1
	obs<-data.frame(cbind(1:n1,rep(0,n1),rep(0,n1),rep(0,n1),rep(0,n1)))
	names(obs)<-c('patient','level','tox','a','order')

#######Beginning at dose level 1
###1st Stage:  Up and down scheme with cohort size 1


i<-1
x0<-lapply(zones,ff)
##'initial.scheme' is a vector indicating the Stage I escalation scheme
initial.scheme<-as.numeric(unlist(x0))
initial.scheme<-c(initial.scheme,rep(d,n-length(initial.scheme)))

while(i < n){
obs$order[i]<-99
obs$level[i]<-initial.scheme[i]
p<-runif(1)
#number of tox in 1st patient
index<-p<=r[obs$level[i]] ##determines any toxicities
if(any(index)){
 obs$tox[i]<-obs$tox[i]+as.numeric(index)
	break
 }###if observe at least 1 toxicity, go to stage 2
i<-i+1
if(length(obs$level[obs$level==d])==stop){
	MTD<-d
	break
}
}
	

##2nd stage
N<<-table(obs$level>0)[2]+1
if(sum(obs$tox[1])==1){
	MTD<-0
	}
	else if(any(obs$tox>0)){
		  level<-crm(obs[1:(N-1),],alpha,theta)
		  obs$a[N-1]<-ahat
	        obs$order[N-1]<-ord

for(j in N:n1){
##assigment for remaining patients
	obs$level[j]<-level
	if(obs$level[n1]>0){
			MTD<-obs$level[n1]
		break
		}

	if(length(obs$level[obs$level==level])==stop+1){
		MTD<-level
		break
		}

	index<-runif(1)<=r[obs$level[j]]
	if(index){obs$tox[j]<-1}
	level<-crm(obs[1:j,],alpha,theta)
	obs$a[j]<-ahat
	obs$order[j]<-ord
	
##crm dose recommendation for Nth patient
	}
	} else
	MTD<-d
	out<-list(trial=obs[obs$level>0,],MTD.rec=MTD)
	}
###'twostgcrm' ENDS HERE


##LOAD FUNCTION 'pocrm.sim'

pocrm.sim<-function(B,accept){
dlt<-p<-ca<-rep(0,d)
rec0<-ss<-rep(0,B)
a<-b<-matrix(,nrow=d,ncol=B)
for(g in 1:B){
	fit<-twostgcrm(n,alpha,r,theta)
	ss[g]<-sum(fit$trial$order>0)
	   for(w in 1:d){
			a[w,g]<-sum(fit$trial$level[1:ss[g]]==w)
			b[w,g]<-sum(fit$trial$level[1:ss[g]]==w & fit$trial$tox[1:ss[g]]==1)		
			rec0[g]<-fit$MTD.rec
			}
}
for(h in 1:d){
	p[h]<-sum(rec0==h)/B
	ca[h]<-sum(a[h,])/sum(ss)
}
output<-list(true.prob=r,MTD.selection=round(p,2),patient.allocation=round(ca,2),percent.DLT=sum(b)/sum(ss),mean.n=mean(ss[ss>1]),stopped=1-sum(p),acceptable=sum(p[which(round(abs(r-theta),2)<=accept)]))
}
##'pocrm.sim' END HERE


#######################################################
#									
#  This program labels the dose- schedule combinations
#  according to the rows of the matrix.
#  For instance, for a 4 x 3 matrix, the labels 
#  would be...
#
#	10	11	12
#	 7	 8	 9
#	 4	 5 	 6
#	 1 	 2 	 3
#
#  where combination 1 corresponds to the lowest 
#  dose and lowest schedule, and so on...
#
######################################################


#######################################################
#									
#  Input:
#	
######################################################

#####number of dose-schedule combinations
d<-12

#####number of possible orderings
s<-12

###Specifiy the possible orderings of the combinations
orders<-matrix(nrow=s,ncol=d)
orders[1,]<-c(1,2,3,4,5,6,7,8,9,10,11,12)
orders[2,]<-c(1,2,4,3,5,7,6,8,10,9,11,12)
orders[3,]<-c(1,4,2,7,5,3,10,8,6,11,9,12)
orders[4,]<-c(1,4,7,10,2,5,8,11,3,6,9,12)
orders[5,]<-c(1,2,4,7,5,3,6,8,10,11,9,12)
orders[6,]<-c(1,4,2,3,5,7,10,8,6,9,11,12)

orders[7,]<-c(1,2,3,7,8,9,4,5,6,10,11,12)
orders[8,]<-c(1,7,4,10,2,8,5,11,3,9,6,12)
orders[9,]<-c(1,2,7,3,8,4,9,5,10,6,11,12)
orders[10,]<-c(1,7,2,4,8,3,10,5,9,11,6,12)
orders[11,]<-c(1,2,7,4,8,3,9,5,10,11,6,12)
orders[12,]<-c(1,7,2,3,8,4,10,5,9,6,11,12)

###specify the zones for use in Stage I
zones<-list(z1=c(1),z2=c(2,4,7),z3=c(3,8,5),z4=c(6,10,9),z5=c(11),z6=c(d))
ff<-function(x){
	if(length(x)==1){
		x
	} else
	sample(x)
}

####specify the target DLT rate
theta<-0.30

####specify the maximum sample size
n<-60

###specify the number consecutive patients on 
###a combination needed to stop the trial
stop<-61

###Specify 'initial guesses' of DLT probabilities (called 'skeleton') for each combination 
###User can use the function 'getprior' in the package 'dfcrm'
###'getprior' implements the algorithm of Lee SM, Cheung, YK.Model calibration in the continual reassessment method. Clin
###Trials 2009 Jun;6(3):227-38.
###'Reasonable' initial guesses will have adequate spacing between values at adjacent levels.
###For more information see O’Quigley J, Zohar, S. Retrospective robustness of the continual reassessment method.
###Journal of Biopharmaceutical Statistics 2010;20(5):1013-25.

library(dfcrm)
skeleton<-round(getprior(0.05,theta,6,d),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
alpha<-matrix(0,nrow=s,ncol=d)
for(j in 1:s){
	alpha[j,]<-skeleton[order(orders[j,])]
}

#######################################################
#									
#  Simulate a single trial using 'twostgcrm'
#	
######################################################

###'trial' returns:
###		'patient' = patient number
###		'level' = combination given to each patient
###		'tox' = DLT indicator for each patient
###		'a' = parameter estimate of a
###		'order' = estimate of the toxicity ordering

###'MTD.rec' returns the combination estimated to be the MTD

fit<-twostgcrm(n,alpha,r,theta)
fit


#######################################################
#									
#  Simulate 'B' trials using 'pocrm.sim'
#	
######################################################

###'fit.sim' returns

###'true.prob' is the true DLT rates
###'MTD.selection' is the proportion of times each combo is recommended as MTD
###'patient.allocation' is the proportion of patients allocated to each combination
###'mean.n' is the average trial size
###'stopped' is the proportion of trials stopped early for safety
###'acceptable' is the proportion of MTD recommendation within +/- 'accept' percent of the target

r1<-c(0.05,0.07,0.11,0.09,0.12,0.18,0.16,0.18,0.23,0.22,0.26,0.30)
r2<-c(0.03,0.14,0.28,0.09,0.21,0.40,0.18,0.32,0.54,0.31,0.45,0.62)
r3<-c(0.10,0.26,0.35,0.30,0.32,0.50,0.45,0.55,0.62,0.60,0.62,0.72)
r4<-c(0.50,0.54,0.58,0.53,0.60,0.65,0.55,0.65,0.75,0.57,0.73,0.78)
r5<-c(0.10,0.28,0.45,0.12,0.30,0.48,0.14,0.32,0.55,0.30,0.48,0.70)
r6<-c(0.03,0.15,0.30,0.12,0.30,0.50,0.30,0.50,0.60,0.50,0.60,0.75)
r7<-c(0.01,0.10,0.50,0.03,0.30,0.55,0.05,0.50,0.60,0.10,0.60,0.70)

r<-r6
fit.sim<-pocrm.sim(B=1000,accept=0.10)
fit.sim