#######################################################
#
#
#   R code for simulating the optimal design for 
#   estimating the most successful dose (MSD)
#   Zohar and O'Quigley (Stats in Med, 2006).
#
#
#######################################################



msdoptimal<-function(r,q,n,B){
select=rep(0,length(r))
j=1
while(j<=B){

	latent.tox=runif(n,0,1)
	latent.eff=runif(n,0,1)

	tox=eff=matrix(,nrow=n,ncol=length(r))
	i=1
	while(i<=n){
		tox[i,]<-as.numeric(latent.tox[i]<=r)
		eff[i,]<-as.numeric(latent.eff[i]<=q)
		i=i+1
	}	
	utility=colMeans(eff)*(1-colMeans(tox))
	sugglev=which(utility==max(utility))
	msd=ifelse(length(sugglev)==1,sugglev,sample(sugglev,1))
	select[msd]=select[msd]+1
	j=j+1
}
cat("True tox probability: ",     round(r,3), sep="    ",  "\n");
cat("True eff probability: ",     round(q,3), sep="    ",  "\n");
cat("True suc probability: ", round(q*(1-r),2), sep="    ",  "\n");
cat("Selection percentage: ", formatC((select/B)*100, digits=1, format="f"), sep="    ",  "\n");
}


##True scenarios
r1=c(0.05,0.15,0.30,0.40,0.50,0.60)
q1=c(0.28,0.36,0.45,0.52,0.65,0.79)

r2=c(0.30,0.38,0.45,0.55,0.70,0.80)
q2=c(0.45,0.55,0.65,0.75,0.85,0.95)

r3=c(0.18,0.28,0.36,0.44,0.52,0.65)
q3=c(0.40,0.50,0.60,0.70,0.80,0.90)


####################################
#
#	Input:
#
#	r = true toxicity probabilities 
#	n = sample size
#	q = true efficacy probabilities 
#	B = number of simulated trials
# 
####################################
n=25
B=10000
r=r1
q=q1
msdoptimal(r,q,n,B)