#rm(list=ls())
# setwd("C:/Users/gzhou/Dropbox/Academic/RA/SeqDesign")
# setwd("/home/guohai/Dropbox/Academic/RA/SeqDesign/")
library(Iso)

intercept = log(0.9/0.1)-21.01842*0.85
p_curve = data.frame(coef=c(intercept  ,  21.01842 ))
# log(0.9/0.1) = x + 21.01842*0.85

curve_eq = function(x) 1/( 1 + exp(-p_curve$coef[1]-p_curve$coef[2]*x) )
true_ed90 = as.numeric( ( log(0.9/0.1)-p_curve$coef[1])/p_curve$coef[2] )
true_ed50 = as.numeric( ( log(0.5/0.5)-p_curve$coef[1])/p_curve$coef[2] )
curve(curve_eq,from=0.7,to=1.1) # true dose-response curve

all_dose_seq = seq(0.7,1,by=0.05)
 
# a1+(4-1)*0.05 = 0.850
# 0.7+(7-1)*0.05  
true_prior_dose_seq = all_dose_seq[2:6]
true_prior = curve_eq(true_prior_dose_seq)

non_informative_prior_dose_seq = c(0.4,0.55,0.70,0.85,1.0)
non_informative_prior_line = function(x) { slope = (true_prior[5]-0.90)/(true_prior_dose_seq[5]-0.85);  slope*(x-0.85)+0.90   }
non_informative_prior = c(0.88, 0.89, 0.90, 0.91, 0.92)

high_prior_dose_seq = all_dose_seq[1:5]
high_prior = curve_eq(high_prior_dose_seq)*0.900/curve_eq(high_prior_dose_seq)[3]
high_prior[high_prior>1] = 0.99

low_prior_dose_seq = all_dose_seq[3:7]
low_prior = curve_eq(low_prior_dose_seq)*0.900/curve_eq(low_prior_dose_seq)[3]



library(bcrm);
set.seed(12);
#nsim = 1000
nsim = 60
target_prob = 0.9
prob_at_ED90_stop_CI = c(0.82,0.99)
prior.alpha = list(1,1,1)
max_trial_size = 60


#### the specification of this dose sequence reflects some clinical belief: the true ED90 is somewhere between 0.7 and 0.9, consequently the interval from 0.7 to 0.9 was more finely split. 

#### Since CRM ALWAYS concludes one of the specified dose levels to be the final ED90, the specification of the dose sequence will be important.  If the dose sequence fails to cover the true ED90 (which is something clinicians must prevent), the performance will be inevitably poor.


#### I've tried to match the CRM and the BCD simulation settings: a) the common maximum sample size is 60
# b) a common true dose-curve for data-generation is used

###### visualize prior
# pdf("figure_4_col_June2016.pdf",width=9,height=6)
range_prob = range(c(true_prior,non_informative_prior,low_prior,high_prior ))
setEPS(); postscript("no_box_figure4_21oct16.eps",width=9,height=6)

plot(true_prior~true_prior_dose_seq,xlim=c(0.4,1),ylim=c(0.30,1),xlab="etHAL dose level [%]",ylab="Probability of response",type="o",lwd=3,pch=0,lty=1,cex=1.5,col="black",axes=F) #  # xlab="HAL dose"
axis(1,seq(0.4,1,by=0.1) )
axis(2,seq(0.3,1,by=0.1))
abline(h=0.90,lty=6,lwd=0.6)
points(low_prior_dose_seq,low_prior,pch=1,cex=1.5,col="red")
lines(low_prior_dose_seq,low_prior,lty=2,lwd=3,col="red")
points(high_prior_dose_seq,high_prior,pch=2,cex=1.5,col="blue")
lines(high_prior_dose_seq,high_prior,lty=3,lwd=3,col="blue")
points(non_informative_prior_dose_seq,non_informative_prior,pch=5,cex=1.5,col="orange")
lines(non_informative_prior_dose_seq,non_informative_prior,lty=4,lwd=3,col="orange")
legend("bottomright",c("True","Low","High","Non-info"),lty=c(1,2,3,4),pch=c(0,1,2,5),lwd=2,cex=1.7,col=c("black","red","blue","orange"))



cohort.size = 2; 
power_m2_low_prior = bcrm(stop = list(nmax = 0, precision =  prob_at_ED90_stop_CI), p.tox0 = low_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= low_prior_dose_seq, nsims = 1, truep = curve_eq(low_prior_dose_seq), quietly = FALSE)
power_m2_low_prior[[1]]$ndose[[1]]$quantiles[c(1,5),] # 95% of CI based on prior
power_m2_low_prior[[1]]$ndose[[1]]$est # estimate with 0 sample size are just prior guess
dev.off()


## a function to extract sample size and RMSE, for a given simulation outcome object (that is returned by bcrm)
simulation_summary = function(bcrm_outcome_obj, nsim1=nsim)
{
  estimated_prob_at_ED90 = rep(NA,nsim1);  estimated_ED90 = rep(NA,nsim1);
  sample_size = rep(NA,nsim1)
  for(k in 1:nsim1)
  { 
    #bcrm_outcome_obj[[k]] is a list including the simulation details for the k-th simulated trial 
    
    
    if(  !( "data"%in%names(bcrm_outcome_obj[[k]]) )  ) #special case, prior guess already has 95% CI within the required precision, CRM converges at the beginning, in this special case, the matrix bcrm_outcome_obj[[k]]$data does NOT exist
    { sample_size[k] = 0; }
    
    else {
      sample_size[k] = nrow(bcrm_outcome_obj[[k]]$data) #  
    }
    ncohort = sample_size[k]/cohort.size
    
    
    final_est_dose_index = bcrm_outcome_obj[[k]]$ndose[[ncohort+1]]$ndose 
    # the index for the dose that was considered as the final ED90 estimate; 
    
    estimated_ED90[k] = bcrm_outcome_obj[[k]]$dose[final_est_dose_index] 
    
    estimated_prob_at_ED90[k] = bcrm_outcome_obj[[k]]$ndose[[ncohort+1]]$est[final_est_dose_index]  
    
  }
  sample_size = sample_size
  rmse_dose = sqrt( mean( ( (estimated_ED90-true_ed90)^2 ) ) )
  rmse_prob_with_target_prob = sqrt( mean( ( (estimated_prob_at_ED90-target_prob)^2 ) ) )
  rmse_prob = sqrt( mean( ( (estimated_prob_at_ED90-  curve_eq(estimated_ED90)   )^2 ) ) )  # root_mean_squared_prob is defined relative to true_probabilty at estimated_ED90 rather than 90%; only when estimated_ED90 = true_ED90 do we have  curve_eq(estimated_ED90) =  target_prob =90%
  list(sample_size = sample_size, rmse_dose=rmse_dose, rmse_prob=rmse_prob,rmse_prob_with_target_prob=rmse_prob_with_target_prob,estimated_ED90=estimated_ED90,estimated_prob_at_ED90=estimated_prob_at_ED90)
}



####*********** cohort size 2, working model power, low_prior 
cohort.size = 2; 
simu_object_power_m2_low_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = low_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= low_prior_dose_seq, nsims = nsim, truep = curve_eq(low_prior_dose_seq), quietly = FALSE)
summary_power_m2_low_prior = simulation_summary(simu_object_power_m2_low_prior)




####*********** cohort size 2, working model logistic, low_prior 
cohort.size = 2; 
simu_object_logistic_m2_low_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = low_prior, ff="logit1",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= low_prior_dose_seq, nsims = nsim, truep = curve_eq(low_prior_dose_seq), quietly = FALSE)
summary_logistic_m2_low_prior = simulation_summary(simu_object_logistic_m2_low_prior)



####*********** cohort size 2, working model power, high_prior 
cohort.size = 2;  
simu_object_power_m2_high_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = high_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose = high_prior_dose_seq, nsims = nsim, truep = curve_eq(high_prior_dose_seq), quietly = FALSE)
summary_power_m2_high_prior = simulation_summary(simu_object_power_m2_high_prior)



####*********** cohort size 2, logistic model power, high_prior 
cohort.size = 2;  
simu_object_logistic_m2_high_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = high_prior, ff="logit1",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= high_prior_dose_seq, nsims = nsim, truep = curve_eq(high_prior_dose_seq), quietly = FALSE)
summary_logistic_m2_high_prior = simulation_summary(simu_object_logistic_m2_high_prior)


####*********** cohort size 2, working model power, true_prior 
cohort.size = 2; 
simu_object_power_m2_true_prior = bcrm(stop = list(nmax = max_trial_size, precision = prob_at_ED90_stop_CI), p.tox0 = true_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= true_prior_dose_seq, nsims = nsim, truep = curve_eq(true_prior_dose_seq), quietly = FALSE)
summary_power_m2_true_prior = simulation_summary(simu_object_power_m2_true_prior)



####*********** cohort size 2, working model logistic, true_prior 
cohort.size = 2; 
simu_object_logistic_m2_true_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = true_prior, ff="logit1",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= true_prior_dose_seq, nsims = nsim, truep = curve_eq(true_prior_dose_seq), quietly = FALSE)
summary_logistic_m2_true_prior = simulation_summary(simu_object_logistic_m2_true_prior)



####*********** cohort size 2, working model power, non_informative_prior 
cohort.size = 2; 
simu_object_power_m2_non_informative_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = non_informative_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= non_informative_prior_dose_seq, nsims = nsim, truep = curve_eq(non_informative_prior_dose_seq), quietly = FALSE)
summary_power_m2_non_informative_prior = simulation_summary(simu_object_power_m2_non_informative_prior)



####*********** cohort size 2, working model logistic, non_informative_prior 
cohort.size = 2; 
simu_object_logistic_m2_non_informative_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = non_informative_prior, ff="logit1",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= non_informative_prior_dose_seq, nsims = nsim, truep = curve_eq(non_informative_prior_dose_seq), quietly = FALSE)
summary_logistic_m2_non_informative_prior = simulation_summary(simu_object_logistic_m2_non_informative_prior)























####*********** cohort size 3, working model power, low_prior 
cohort.size = 3; 
simu_object_power_m3_low_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = low_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= low_prior_dose_seq, nsims = nsim, truep = curve_eq(low_prior_dose_seq), quietly = FALSE)
summary_power_m3_low_prior = simulation_summary(simu_object_power_m3_low_prior)




####*********** cohort size 3, working model logistic, low_prior 
cohort.size = 3; 
simu_object_logistic_m3_low_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = low_prior, ff="logit1",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= low_prior_dose_seq, nsims = nsim, truep = curve_eq(low_prior_dose_seq), quietly = FALSE)
summary_logistic_m3_low_prior = simulation_summary(simu_object_logistic_m3_low_prior)



####*********** cohort size 3, working model power, high_prior 
cohort.size = 3;  
simu_object_power_m3_high_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = high_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= high_prior_dose_seq, nsims = nsim, truep = curve_eq(high_prior_dose_seq), quietly = FALSE)
summary_power_m3_high_prior = simulation_summary(simu_object_power_m3_high_prior)



####*********** cohort size 3, working model logistic, high_prior
cohort.size = 3;  
simu_object_logistic_m3_high_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = high_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= high_prior_dose_seq, nsims = nsim, truep = curve_eq(high_prior_dose_seq), quietly = FALSE)
summary_logistic_m3_high_prior = simulation_summary(simu_object_logistic_m3_high_prior)


####*********** cohort size 3, working model power, true_prior 
cohort.size = 3; 
simu_object_power_m3_true_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = true_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= true_prior_dose_seq, nsims = nsim, truep = curve_eq(true_prior_dose_seq), quietly = FALSE)
summary_power_m3_true_prior = simulation_summary(simu_object_power_m3_high_prior)



####*********** cohort size 3, working model logistic, true_prior
cohort.size = 3; 
simu_object_logistic_m3_true_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = true_prior, ff="logit1",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= true_prior_dose_seq, nsims = nsim, truep = curve_eq(true_prior_dose_seq), quietly = FALSE)
summary_logistic_m3_true_prior = simulation_summary(simu_object_logistic_m3_true_prior)



####*********** cohort size 3, working model power, non_informative_prior 
cohort.size = 3; 
simu_object_power_m3_non_informative_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = non_informative_prior, ff="power",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= non_informative_prior_dose_seq, nsims = nsim, truep = curve_eq(non_informative_prior_dose_seq), quietly = FALSE)
summary_power_m3_non_informative_prior = simulation_summary(simu_object_power_m3_non_informative_prior)



####*********** cohort size 3, working model logistic, non_informative_prior 
cohort.size = 3; 
simu_object_logistic_m3_non_informative_prior = bcrm(stop = list(nmax = max_trial_size, precision =  prob_at_ED90_stop_CI), p.tox0 = non_informative_prior, ff="logit1",cohort = cohort.size, target.tox=target_prob, constrain = FALSE, prior.alpha = prior.alpha, simulate = TRUE,dose= non_informative_prior_dose_seq, nsims = nsim, truep = curve_eq(non_informative_prior_dose_seq), quietly = FALSE)
summary_logistic_m3_non_informative_prior = simulation_summary(simu_object_logistic_m3_non_informative_prior)


save.image("CRM_2016_Oct21")
 
#### create boxplots to summarize simulation results 
# pdf("figure_5_June2016.pdf",width = 10, height=7.3)
setEPS(); postscript("no_box_figure5_21oct2016.eps",width=10,height=7.3)
par(mfrow=c(2,2))
sz = c(summary_power_m2_true_prior$sample_size, summary_power_m2_low_prior$sample_size,  summary_power_m2_high_prior$sample_size, summary_power_m2_non_informative_prior$sample_size) 
grp = c(rep(1,nsim),rep(2,nsim),rep(3,nsim),rep(4,nsim))  
grp = as.factor(grp)
boxplot(sz~grp,main="Power working model, cohort size 2",ylab="Sample size",xlab="Type of prior guess of response probability",ylim=c(-10,60),xaxt="n",yaxt="n",axes=F)
axis(1, at = 1:4, labels = c("True prior","Low prior","High prior","Non-info prior"))
axis(2,at=seq(-10,60,by=10),labels=c("","","10","20","30","40","50","60"))
for(k in 1:4)
{  
  if(k==1) { RMSE_dose = summary_power_m2_true_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m2_true_prior$rmse_prob; }
  if(k==2) { RMSE_dose = summary_power_m2_low_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m2_low_prior$rmse_prob; }
  if(k==3) { RMSE_dose = summary_power_m2_high_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m2_high_prior$rmse_prob; }
  if(k==4) { RMSE_dose = summary_power_m2_non_informative_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m2_non_informative_prior$rmse_prob; }  
  
  text(k,0-4,paste("RMSE_d=",as.character(round(RMSE_dose,2)),sep=""),cex=0.77)
  text(k,0-8,paste("RMSE_p=",as.character(round(RMSE_dose_prob_at_ED90,2)),sep=""),cex=0.77)
}


sz = c(summary_power_m3_true_prior$sample_size, summary_power_m3_low_prior$sample_size,  summary_power_m3_high_prior$sample_size, summary_power_m3_non_informative_prior$sample_size) 
grp = c(rep(1,nsim),rep(2,nsim),rep(3,nsim),rep(4,nsim))  
grp = as.factor(grp)
boxplot(sz~grp,main="Power working model, cohort size 3",ylab="Sample size",xlab="Type of prior guess of response probability", ylim=c(-10,max_trial_size),xaxt="n",yaxt="n",axes=F)
axis(1, at = 1:4, labels = c("True prior","Low prior","High prior","Non-info prior"))
axis(2,at=seq(-10,60,by=10),labels=c("","","10","20","30","40","50","60"))
for(k in 1:4)
{  
  if(k==1) { RMSE_dose = summary_power_m3_true_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m3_true_prior$rmse_prob; }
  if(k==2) { RMSE_dose = summary_power_m3_low_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m3_low_prior$rmse_prob; }
  if(k==3) { RMSE_dose = summary_power_m3_high_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m3_high_prior$rmse_prob; }
  if(k==4) { RMSE_dose = summary_power_m3_non_informative_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_power_m3_non_informative_prior$rmse_prob; }  
  
  text(k,0-4,paste("RMSE_d=",as.character(round(RMSE_dose,2)),sep=""),cex=0.77)
  text(k,0-8,paste("RMSE_p=",as.character(round(RMSE_dose_prob_at_ED90,2)),sep=""),cex=0.77)
}




sz = c(summary_logistic_m2_true_prior$sample_size, summary_logistic_m2_low_prior$sample_size,  summary_logistic_m2_high_prior$sample_size, summary_logistic_m2_non_informative_prior$sample_size) 
grp = c(rep(1,nsim),rep(2,nsim),rep(3,nsim),rep(4,nsim))  
grp = as.factor(grp)
boxplot(sz~grp,main="Logistic working model, cohort size 2",ylab="Sample size",xlab="Type of prior guess of response probability",
        ylim=c(-10,max_trial_size),xaxt="n",yaxt="n",axes=F)
axis(1, at = 1:4, labels = c("True prior","Low prior","High prior","Non-info prior"))
axis(2,at=seq(-10,60,by=10),labels=c("","","10","20","30","40","50","60"))
for(k in 1:4)
{  
  if(k==1) { RMSE_dose = summary_logistic_m2_true_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m2_true_prior$rmse_prob; }
  if(k==2) { RMSE_dose = summary_logistic_m2_low_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m2_low_prior$rmse_prob; }
  if(k==3) { RMSE_dose = summary_logistic_m2_high_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m2_high_prior$rmse_prob; }
  if(k==4) { RMSE_dose = summary_logistic_m2_non_informative_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m2_non_informative_prior$rmse_prob; }  
  
  text(k,0-4,paste("RMSE_d=",as.character(round(RMSE_dose,2)),sep=""),cex=0.77)
  text(k,0-8,paste("RMSE_p=",as.character(round(RMSE_dose_prob_at_ED90,2)),sep=""),cex=0.77)
}



sz = c(summary_logistic_m3_true_prior$sample_size, summary_logistic_m3_low_prior$sample_size,  summary_logistic_m3_high_prior$sample_size, summary_logistic_m3_non_informative_prior$sample_size) 
grp = c(rep(1,nsim),rep(2,nsim),rep(3,nsim),rep(4,nsim))  
grp = as.factor(grp)
boxplot(sz~grp,main="Logistic working model, cohort size 3",ylab="Sample size",xlab="Type of prior guess of response probability", ylim=c(-10,max_trial_size),xaxt="n",yaxt="n",axes=F)
axis(1, at = 1:4, labels = c("True prior","Low prior","High prior","Non-info prior"))
axis(2,at=seq(-10,60,by=10),labels=c("","","10","20","30","40","50","60"))
for(k in 1:4)
{  
  if(k==1) { RMSE_dose = summary_logistic_m3_true_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m3_true_prior$rmse_prob; }
  if(k==2) { RMSE_dose = summary_logistic_m3_low_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m3_low_prior$rmse_prob; }
  if(k==3) { RMSE_dose = summary_logistic_m3_high_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m3_high_prior$rmse_prob; }
  if(k==4) { RMSE_dose = summary_logistic_m3_non_informative_prior$rmse_dose; RMSE_dose_prob_at_ED90 = summary_logistic_m3_non_informative_prior$rmse_prob; }  
  
  text(k,0-4,paste("RMSE_d=",as.character(round(RMSE_dose,2)),sep=""),cex=0.77)
  text(k,0-8,paste("RMSE_p=",as.character(round(RMSE_dose_prob_at_ED90,2)),sep=""),cex=0.77)
}
dev.off()

rm(list=ls())
load("CRM_2016_Oct21")
## export 3 matrices for Rollin's code to plot: matrix of sample size, matrix of rmse_prob, matrix of rmse_dose
sample_size_matrix = matrix( c(summary_power_m2_true_prior$sample_size,summary_power_m2_low_prior$sample_size,summary_power_m2_high_prior$sample_size,summary_power_m2_non_informative_prior$sample_size, summary_power_m3_true_prior$sample_size,summary_power_m3_low_prior$sample_size,summary_power_m3_high_prior$sample_size,summary_power_m3_non_informative_prior$sample_size,  summary_logistic_m2_true_prior$sample_size,summary_logistic_m2_low_prior$sample_size,summary_logistic_m2_high_prior$sample_size,summary_logistic_m2_non_informative_prior$sample_size, summary_logistic_m3_true_prior$sample_size,summary_logistic_m3_low_prior$sample_size,summary_logistic_m3_high_prior$sample_size,summary_logistic_m3_non_informative_prior$sample_size),ncol=16,dimnames = list(paste0("simu_dataset_No.",1:nsim), c("True_prior_power_model_cohort2","Low_prior_power_model_cohort2","High_prior_power_model_cohort2","non_info_prior_power_model_cohort2","True_prior_power_model_cohort3","Low_prior_power_model_cohort3","High_prior_power_model_cohort3","non_info_prior_power_model_cohort3","True_prior_logistic_model_cohort2","Low_prior_logistic_model_cohort2","High_prior_logistic_model_cohort2","non_info_prior_logistic_model_cohort2","True_prior_logistic_model_cohort3","Low_prior_logistic_model_cohort3","High_prior_logistic_model_cohort3","non_info_logistic_power_model_cohort3")))


rmse_prob_matrix = matrix(c(summary_power_m2_true_prior$rmse_prob,summary_power_m2_low_prior$rmse_prob,summary_power_m2_high_prior$rmse_prob,summary_power_m2_non_informative_prior$rmse_prob, summary_power_m3_true_prior$rmse_prob,summary_power_m3_low_prior$rmse_prob,summary_power_m3_high_prior$rmse_prob,summary_power_m3_non_informative_prior$rmse_prob, summary_logistic_m2_true_prior$rmse_prob,summary_logistic_m2_low_prior$rmse_prob,summary_logistic_m2_high_prior$rmse_prob,summary_logistic_m2_non_informative_prior$rmse_prob, summary_logistic_m3_true_prior$rmse_prob,summary_logistic_m3_low_prior$rmse_prob,summary_logistic_m3_high_prior$rmse_prob,summary_logistic_m3_non_informative_prior$rmse_prob),ncol=16,dimnames = list(c("rmse_prob_m"), c("True_prior_power_model_cohort2","Low_prior_power_model_cohort2","High_prior_power_model_cohort2","non_info_prior_power_model_cohort2","True_prior_power_model_cohort3","Low_prior_power_model_cohort3","High_prior_power_model_cohort3","non_info_prior_power_model_cohort3","True_prior_logistic_model_cohort2","Low_prior_logistic_model_cohort2","High_prior_logistic_model_cohort2","non_info_prior_logistic_model_cohort2","True_prior_logistic_model_cohort3","Low_prior_logistic_model_cohort3","High_prior_logistic_model_cohort3","non_info_logistic_power_model_cohort3")))

rmse_dose_matrix = matrix(c(summary_power_m2_true_prior$rmse_dose,summary_power_m2_low_prior$rmse_dose,summary_power_m2_high_prior$rmse_dose,summary_power_m2_non_informative_prior$rmse_dose, summary_power_m3_true_prior$rmse_dose,summary_power_m3_low_prior$rmse_dose,summary_power_m3_high_prior$rmse_dose,summary_power_m3_non_informative_prior$rmse_dose, summary_logistic_m2_true_prior$rmse_dose,summary_logistic_m2_low_prior$rmse_dose,summary_logistic_m2_high_prior$rmse_dose,summary_logistic_m2_non_informative_prior$rmse_dose, summary_logistic_m3_true_prior$rmse_dose,summary_logistic_m3_low_prior$rmse_dose,summary_logistic_m3_high_prior$rmse_dose,summary_logistic_m3_non_informative_prior$rmse_dose),ncol=16,dimnames = list(c("rmse_prob_m"), c("True_prior_power_model_cohort2","Low_prior_power_model_cohort2","High_prior_power_model_cohort2","non_info_prior_power_model_cohort2","True_prior_power_model_cohort3","Low_prior_power_model_cohort3","High_prior_power_model_cohort3","non_info_prior_power_model_cohort3","True_prior_logistic_model_cohort2","Low_prior_logistic_model_cohort2","High_prior_logistic_model_cohort2","non_info_prior_logistic_model_cohort2","True_prior_logistic_model_cohort3","Low_prior_logistic_model_cohort3","High_prior_logistic_model_cohort3","non_info_logistic_power_model_cohort3")))

save(list=c("sample_size_matrix","rmse_prob_matrix","rmse_dose_matrix"),file="crm_simul.rda")