####################################################### # # # 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)