diff --git a/DESCRIPTION b/DESCRIPTION
index 8be0ed4..e15756a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,26 +2,29 @@ Package: BayesianMCPMod
 Title: Bayesian MCPMod
 Version: 0.1.3-3
 Authors@R: c(
-    person("Sebastian", "Bossert",
-           role = "aut",
-           email = "sebastian.bossert@boehringer-ingelheim.com"),
+    person("Boehringer Ingelheim Pharma GmbH & Co. KG",
+           role = c("cph", "fnd")),
     person("Stephan", "Wojciekowski",
            role = c("aut", "cre"),
            email = "stephan.wojciekowski@boehringer-ingelheim.com"),
     person("Lars", "Andersen",
            role = "aut",
            email = "lars.andersen@boehringer-ingelheim.com"),
-    person("Boehringer Ingelheim Pharma GmbH & Co. KG",
-           role = c("cph", "fnd"))
+    person("Steven", "Brooks",
+           role = "aut",
+           email = "lars.andersen@boehringer-ingelheim.com"),
+    person("Sebastian", "Bossert",
+           role = "aut",
+           email = "sebastian.bossert@boehringer-ingelheim.com")
     )
-Description: Provides functions for the analysis of dose finding trials with Bayesian MCPMod.
+Description: Provides functions for the planning and analysis of dose finding trials with Bayesian MCPMod.
 License: Apache License (>= 2)
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Imports:
-    DoseFinding,
+    DoseFinding (>= 1.1-1),
     ggplot2,
     stats,
     RBesT,
diff --git a/NAMESPACE b/NAMESPACE
index f89ff4d..dc6a5e4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,19 +1,19 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(plot,modelFits)
-S3method(predict,ModelFits)
+S3method(predict,modelFits)
 S3method(print,BayesianMCP)
 S3method(print,BayesianMCPMod)
 S3method(print,modelFits)
 S3method(print,postList)
 S3method(summary,postList)
 export(assessDesign)
-export(getBootstrapBands)
-export(getContrMat)
+export(getBootstrapQuantiles)
+export(getContr)
 export(getCritProb)
+export(getESS)
 export(getModelFits)
 export(getPosterior)
-export(getPriorList)
 export(performBayesianMCP)
 export(performBayesianMCPMod)
 export(simulateData)
diff --git a/R/BMCPMod.R b/R/BMCPMod.R
index 66522e0..4dafc4e 100644
--- a/R/BMCPMod.R
+++ b/R/BMCPMod.R
@@ -3,12 +3,37 @@
 #' @description This function performs simulation based trial design evaluations for a set of specified dose-response models
 #'
 #' @param n_patients Vector specifying the planned number of patients per dose group
-#' @param mods An object of class "Mods" as specified in the Dosefinding package.
+#' @param mods An object of class "Mods" as specified in the DoseFinding package.
 #' @param prior_list a prior_list object specifying the utilized prior for the different dose groups 
+#' @param sd a positive value, specification of assumed sd 
 #' @param n_sim number of simulations to be performed
-#' @param alpha_crit_val critical value to be used for the testing (on the probability scale)
+#' @param alpha_crit_val (unadjusted) critical value to be used for the MCT testing step. Passed to the getCritProb function for the calculation of adjusted critical values (on the probability scale). Default is 0.05.
 #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.
+#' @param reestimate boolean variable, defining whether critical value should be calculated with re-estimated contrasts (see getCritProb function for more details). Default FALSE
+#' @param contr Allows specification of a fixed contrasts matrix. Default NULL
+#' @param dr_means a vector, allows specification of  individual (not model based) assumed effects per dose group. Default NULL
 #' 
+#' @return returns success probabilities for the different assumed dose-response shapes, attributes also includes information around average success rate (across all assumed models) and prior Effective sample size
+#' 
+#' @examples
+#' # example code
+#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6)
+#' sd   = 12
+#' prior_list<-list(Ctrl   = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2),
+#'                    DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+#'                    DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+#'                    DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+#'                    DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+#' 
+#' n_patients <- c(40,60,60,60,60)
+#' success_probabilities <- assessDesign(
+#' n_patients  = n_patients,
+#' mods        = mods,
+#' prior_list  = prior_list,
+#' sd          = sd)
+#' 
+#' success_probabilities
+
 #' @export
 assessDesign <- function (
     
@@ -16,22 +41,45 @@ assessDesign <- function (
   mods,
   prior_list,
   
+  sd,
+  
   n_sim          = 1e3,
   alpha_crit_val = 0.05,
-  simple         = TRUE
+  simple         = TRUE,
+  reestimate     = FALSE,
+  
+  contr          = NULL,
+  dr_means       = NULL
   
 ) {
   
-  dose_levels <- attr(prior_list, "dose_levels")
+  dose_levels <- attr(mods, "doses")
   
   data <- simulateData(
     n_patients  = n_patients,
     dose_levels = dose_levels,
-    sd          = attr(prior_list, "sd_tot"),
+    sd          = sd,
     mods        = mods,
-    n_sim       = n_sim)
+    n_sim       = n_sim, 
+    dr_means    = dr_means)
+  
+  model_names <- colnames(data)[-c(1:3)]
   
-  model_names <- names(mods)
+  crit_prob_adj <- getCritProb(
+    mods           = mods,
+    dose_levels    = dose_levels,
+    dose_weights   = n_patients,
+    alpha_crit_val = alpha_crit_val)
+  
+  if (!reestimate & is.null(contr)) {
+    
+    contr <- getContr(
+      mods           = mods,
+      dose_levels    = dose_levels,
+      dose_weights   = n_patients,
+      prior_list     = prior_list)
+    
+  }
   
   eval_design <- lapply(model_names, function (model_name) {
     
@@ -39,71 +87,173 @@ assessDesign <- function (
       data       = getModelData(data, model_name),
       prior_list = prior_list)
     
-    crit_pval <- getCritProb(
-      mods           = mods,
-      dose_levels    = dose_levels,
-      dose_weights   = n_patients,
-      alpha_crit_val = alpha_crit_val)
-    
-    contr_mat_prior <- getContrMat(
-      mods           = mods,
-      dose_levels    = dose_levels,
-      dose_weights   = n_patients,
-      prior_list     = prior_list)
+    if (reestimate & is.null(contr)) {
+      
+      post_sds <- sapply(posterior_list, function (post) summary(post)[, 2])
+      
+      contr <- apply(post_sds, 2, function (post_sd) getContr(
+        mods           = mods,
+        dose_levels    = dose_levels,
+        sd_posterior   = post_sd))
+      
+    } 
     
     b_mcp_mod <- performBayesianMCPMod(
-      posteriors_list = posterior_list,
-      contr_mat       = contr_mat_prior,
-      crit_prob       = crit_pval,
-      simple          = simple)
+      posterior_list = posterior_list,
+      contr          = contr,
+      crit_prob_adj  = crit_prob_adj,
+      simple         = simple)
     
   })
   
+  avg_success_rate <- mean(sapply(eval_design, function (bmcpmod) {
+    attr(bmcpmod$BayesianMCP, "successRate")
+  }))
+  
   names(eval_design) <- model_names
   
+  attr(eval_design, "avgSuccessRate") <- avg_success_rate
+  attr(eval_design, "placEff")        <- ifelse(test = is.null(dr_means),
+                                                yes  = attr(mods, "placEff"),
+                                                no   = dr_means[1])
+  attr(eval_design, "maxEff")         <- ifelse(test = is.null(dr_means),
+                                                yes  = attr(mods, "maxEff"),
+                                                no   = diff(range(dr_means)))
+  attr(eval_design, "sampleSize")     <- n_patients
+  attr(eval_design, "priorESS")       <- round(getESS(prior_list), 1)
+  
   return (eval_design)
   
 }
 
-#' @title getContrMat
+#' @title getContr
+#' 
+#' @description This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr of the DoseFinding package.
+#' Hereby 4 different options can be distinguished that are automatically executed based on the input that is provided
+#' i)  Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated. 
+#'     In detail,  in a first step the dose_weights (typically the number of patients per dose group) and the prior information is combined by calculating for
+#'     each dose group a posterior effective sample. Based on this posterior effective sample sizes the allocation ratio is derived, which allows for a calculation on
+#'     pseudo-optimal contrasts via regular MCPMod.are calculated from the
+#'     regular MCPMod for these specific weights 
+#' ii) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
+#'     regular MCPMod for these specific weights
+#' iii)Bayesian approach + re-estimation: If only a sd_posterior (i.e. variability of the posterior distribution) is provided, pseudo-optimal contrasts based on these posterior weights will be calculated
+#' iv) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the
+#'     regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances 
 #' 
-#' @description This function calculates contrast vectors that are optimal for detecting certain alternatives. More information and link to publication will be added.
+#' @param mods An object of class "Mods" as specified in the DoseFinding package.
+#' @param dose_levels vector containing the different dosage levels.
+#' @param dose_weights Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option i) is planned) these two inputs should be provided on the same scale (e.g. patient numbers).  Default NULL
+#' @param prior_list a prior_list object, only required as input for Option i). Default NULL
+#' @param sd_posterior a vector of positive values with information about the variability of the posterior distribution, only required for Option iii). Default NULL
+#' @param se_new_trial a vector of positive values with information about the observed variability, only required for Option iv). Default NULL
 #' 
-#' @param mods An object of class "Mods" as specified in the Dosefinding package.
-#' @param dose_levels vector containing the different doseage levels.
-#' @param dose_weights Vector specifying weights for the different doses
-#' @param prior_list a prior_list object
+#' @examples
+#' # example code
+#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6)
+#' dose_levels=c(0, 0.5, 2, 4, 8)
+#' sd_posterior   = c(2.8,3,2.5,3.5,4)
+#' contr_mat<- getContr(
+#' mods           = mods,
+#' dose_levels    = dose_levels,
+#' sd_posterior   = sd_posterior) 
 #' 
-#' @return contr_mat Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package.
+#' @return contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the DoseFinding package.
 #' 
 #' @export
-getContrMat <- function (
-  
+getContr <- function (
+    
   mods,
   dose_levels,
-  dose_weights,
-  prior_list
+  dose_weights = NULL,
+  prior_list   = NULL,
+  sd_posterior = NULL,
+  se_new_trial = NULL
   
 ) {
   
-  ess_prior <- suppressMessages(round(unlist(lapply(prior_list, RBesT::ess))))
+  # frequentist & re-estimation
+  if (!is.null(se_new_trial) & 
+      is.null(dose_weights) & is.null(prior_list) & is.null(sd_posterior)) {
+    
+    w <- NULL
+    S <- diag((se_new_trial)^2)
+    
+    # frequentist & no re-estimation
+  } else if (!is.null(dose_weights) & 
+             is.null(se_new_trial) & is.null(prior_list) & is.null(sd_posterior)) {
+    
+    w <- dose_weights
+    S <- NULL
+    
+    # Bayesian & re-estimation
+  } else if (!is.null(sd_posterior) & 
+             is.null(se_new_trial) & is.null(prior_list) & is.null(dose_weights)) {
+    
+    w <- NULL
+    S <- diag((sd_posterior)^2)
+    
+    # Bayesian & no re-estimation
+  } else if (!is.null(dose_weights) & !is.null(prior_list) & 
+             is.null(se_new_trial) & is.null(sd_posterior)) {
+    
+    w <- dose_weights +
+      suppressMessages(round(unlist(lapply(prior_list, RBesT::ess))))
+    S <- NULL
+    
+  } else {
+    
+    stop (paste("Provided combiations of 'se_new_trial',",
+                "'dose_weights', 'prior_list', 'sd_posterior' not allowed.",
+                "See ?getContr for allowed combinations."))
+    
+  }
   
-  contr_mat <- DoseFinding::optContr(
-    models = mods,
-    doses  = dose_levels,
-    w      = dose_weights + ess_prior)
+  if (is.null(w)) {
+    
+    contr <- DoseFinding::optContr(
+      models = mods,
+      doses  = dose_levels,
+      S      = S)
+    
+  } else {
+    
+    contr <- DoseFinding::optContr(
+      models = mods,
+      doses  = dose_levels,
+      w      = w)
+    
+  }
   
-  return (contr_mat)
+  return (contr)
   
 }
 
 #' @title getCritProb
 #' 
-#' @param mods An object of class "Mods" as specified in the Dosefinding package.
+#' @description This function calculates multiplicity adjusted critical values. The critical values are calculated in such a way that
+#'  when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level). 
+#'  Hereby optimal contrasts of the frequentist MCPMod are applied and two options can be distinguished
+#'  i) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
+#'     regular MCPMod for these specific weights and the corresponding critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package.
+#' ii) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the
+#'     regular MCPMod for this specific vector of standard errors. Here as well the critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package.
+#' 
+#' @param mods An object of class "Mods" as specified in the DoseFinding package.
 #' @param dose_levels vector containing the different dosage levels.
-#' @param dose_weights Vector specifying weights for the different doses
+#' @param dose_weights Vector specifying weights for the different doses, only required for Option i). Default NULL
+#' @param se_new_trial a vector of positive values, only required for Option ii). Default NULL
 #' @param alpha_crit_val significance level. Default set to 0.025.
-#' 
+#'
+#' @examples
+#' # example code
+#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8))
+#' dose_levels=c(0, 0.5, 2, 4, 8)
+#' critVal<- getCritProb(
+#'   mods           = mods,
+#'   dose_weights  =c(50,50,50,50,50), #reflecting the planned sample size
+#'   dose_levels    = dose_levels,
+#'   alpha_crit_val = 0.05) 
 #' @return crit_pval multiplicity adjusted critical value on the probability scale.
 #' 
 #' @export
@@ -111,71 +261,111 @@ getCritProb <- function (
     
   mods,
   dose_levels,
-  dose_weights,
+  dose_weights   = NULL,
+  se_new_trial   = NULL,
   alpha_crit_val = 0.025
   
 ) {
   
-  contr_mat <- DoseFinding::optContr(
-    models = mods,
-    doses  = dose_levels,
-    w      = dose_weights)
+  contr <- getContr(mods           = mods,
+                    dose_levels    = dose_levels ,
+                    dose_weights   = dose_weights,
+                    se_new_trial   = se_new_trial)
   
-  crit_pval <- pnorm(DoseFinding:::critVal(
-    corMat      = contr_mat$corMat,
+  crit_prob <- stats::pnorm(DoseFinding::critVal(
+    corMat      = contr$corMat,
     alpha       = alpha_crit_val,
     df          = 0,
     alternative = "one.sided"))
   
-  return (crit_pval)
+  return (crit_prob)
   
 }
 
 #' @title performBayesianMCPMod
 #' 
-#' @description performs bayesian MCP Test step and modelling.
+#' @description performs bayesian MCP Test step and modelling in a combined fashion. See performBayesianMCP function for MCT Test step and getModelFits for the modelling step
 #' 
-#' @param posteriors_list a getPosterior object
-#' @param contr_mat a getContrMat object, contrast matrix to be used for the testing step.
-#' @param crit_prob a getCritProb object
+#' @param posterior_list a getPosterior object with information about the (mixture) posterior distribution per dose group
+#' @param contr a getContrMat object, contrast matrix to be used for the testing step.
+#' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale).
 #' @param simple boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.
+#' @examples
+#' # example code
+#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8))
+#' dose_levels=c(0, 0.5, 2, 4, 8)
+#' sd_posterior   = c(2.8,3,2.5,3.5,4)
+#' contr_mat<- getContr(
+#' mods           = mods,
+#' dose_levels    = dose_levels,
+#' sd_posterior   = sd_posterior)
+#' critVal<- getCritProb(
+#'   mods           = mods,
+#'   dose_weights  =c(50,50,50,50,50), #reflecting the planned sample size
+#'   dose_levels    = dose_levels,
+#'   alpha_crit_val = 0.05)
+#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
+#'                    DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+#'                    DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+#'                    DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+#'                    DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+#' mu<-c(0,1,1.5,2,2.5)
+#' se<-c(5,4,6,7,8)
+#' posterior_list <- getPosterior(
+#'    prior_list = prior_list,
+#'     mu_hat   = mu,
+#'    se_hat   = se)
+#' performBayesianMCPMod(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal,simple = FALSE)
 #' 
 #' @return bmcpmod test result as well as modelling result.
 #' 
 #' @export
 performBayesianMCPMod <- function (
     
-  posteriors_list,
-  contr_mat,
-  crit_prob,
+  posterior_list,
+  contr,
+  crit_prob_adj,
   simple = FALSE
   
 ) {
   
-  if (class(posteriors_list) == "postList") {
+  if (inherits(posterior_list,  "postList")) {
     
-    posteriors_list <- list(posteriors_list)
+    posterior_list <- list(posterior_list)
     
   }
   
-  b_mcp <- performBayesianMCP(
-    posteriors_list = posteriors_list,
-    contr_mat       = contr_mat,
-    crit_prob       = crit_prob)
+  if (inherits(contr, "optContr")) {
+    
+    model_shapes <- colnames(contr$contMat)
+    dose_levels  <- as.numeric(rownames(contr$contMat))
+    
+  } else if (length(contr) == length(posterior_list)) {
+    
+    model_shapes <- colnames(contr[[1]]$contMat)
+    dose_levels  <- as.numeric(rownames(contr[[1]]$contMat))
+    
+  } else {
+    
+    stop ("Argument 'contr' must be of type 'optContr'")
+    
+  }
   
-  model_shapes <- colnames(contr_mat$contMat)
-  dose_levels  <- as.numeric(rownames(contr_mat$contMat))
+  b_mcp <- performBayesianMCP(
+    posterior_list = posterior_list,
+    contr          = contr,
+    crit_prob_adj  = crit_prob_adj)
   
-  fits_list <- lapply(seq_along(posteriors_list), function (i) {
+  fits_list <- lapply(seq_along(posterior_list), function (i) {
     
     if (b_mcp[i, 1]) {
       
-      sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob")
+      sign_models <- b_mcp[i, -c(1, 2)] > attr(b_mcp, "crit_prob_adj")
       
       model_fits  <- getModelFits(
         models      = model_shapes,
         dose_levels = dose_levels,
-        posterior   = posteriors_list[[i]],
+        posterior   = posterior_list[[i]],
         simple      = simple)
       
       model_fits  <- addSignificance(model_fits, sign_models)
@@ -199,7 +389,7 @@ addSignificance <- function (
     
   model_fits,
   sign_models
-
+  
 ) {
   
   names(sign_models) <- NULL
@@ -216,35 +406,81 @@ addSignificance <- function (
   
 }
 
-#' @title BayesianMCP
+#' @title performBayesianMCP
 #' 
-#' @description performs bayesian MCP Test step.
+#' @description performs bayesian MCP Test step, as described in Fleischer et al. (2022).
+#' Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups).
+#' In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale).
+#' @references Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670. doi:10.1002/pst.2193 
+#' @param posterior_list a getPosterior object with information about the (mixture) posterior distribution per dose group 
+#' @param contr a getContrMat object, contrast matrix to be used for the testing step.
+#' @param crit_prob_adj a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)
 #' 
-#' @param posteriors_list a getPosterior object
-#' @param contr_mat a getContrMat object, contrast matrix to be used for the testing step.
-#' @param crit_prob a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)
+#' @examples
+#' # example code
+#' mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8))
+#' dose_levels=c(0, 0.5, 2, 4, 8)
+#' sd_posterior   = c(2.8,3,2.5,3.5,4)
+#' contr_mat<- getContr(
+#' mods           = mods,
+#' dose_levels    = dose_levels,
+#' sd_posterior   = sd_posterior)
+#' critVal<- getCritProb(
+#'   mods           = mods,
+#'   dose_weights  =c(50,50,50,50,50), #reflecting the planned sample size
+#'   dose_levels    = dose_levels,
+#'   alpha_crit_val = 0.05)
+#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
+#'                    DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+#'                    DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+#'                    DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+#'                    DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+#' mu<-c(0,1,1.5,2,2.5)
+#' se<-c(5,4,6,7,8)
+#' posterior_list <- getPosterior(
+#'    prior_list = prior_list,
+#'     mu_hat   = mu,
+#'    se_hat   = se)
+#' performBayesianMCP(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal)
 #' 
-#' @return b_mcp test result 
+#' @return b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance    
 #' 
 #' @export
 performBayesianMCP <- function(
     
-  posteriors_list,
-  contr_mat,
-  crit_prob
+  posterior_list,
+  contr,
+  crit_prob_adj
   
 ) {
+  if (inherits(posterior_list,  "postList")) {
+    
+    posterior_list <- list(posterior_list)
+    
+  }
   
-  if (class(posteriors_list) == "postList") {
+  if (inherits(contr, "optContr")) {
+    
+    b_mcp <- t(sapply(posterior_list, BayesMCPi, contr, crit_prob_adj))
     
-    posteriors_list <- list(posteriors_list)
+  } else {
+    
+    b_mcp <- t(mapply(BayesMCPi, posterior_list, contr, crit_prob_adj))
     
   }
   
-  b_mcp <- t(sapply(posteriors_list, BayesMCPi, contr_mat, crit_prob))
+  class(b_mcp)                 <- "BayesianMCP"
+  attr(b_mcp, "crit_prob_adj") <- crit_prob_adj
+  attr(b_mcp, "successRate")   <- mean(b_mcp[, 1])
+  attr(b_mcp, "ess_avg")       <- ifelse(
+    test = is.na(attr(posterior_list[[1]], "ess")),
+    yes  = numeric(0),
+    no   = rowMeans(sapply(posterior_list, function (posteriors) {
+      
+      attr(posteriors, "ess")
+      
+    })))
   
-  attr(b_mcp, "crit_prob") <- crit_prob
-  class(b_mcp)             <- "BayesianMCP"
   
   return (b_mcp)
   
@@ -253,8 +489,8 @@ performBayesianMCP <- function(
 BayesMCPi <- function (
     
   posterior_i,
-  contr_mat,
-  crit_prob
+  contr,
+  crit_prob_adj
   
 ) {
   
@@ -282,12 +518,13 @@ BayesMCPi <- function (
   }
   
   post_combs_i <- getPostCombsI(posterior_i)
-  post_probs   <- apply(contr_mat$contMat, 2, getPostProb, post_combs_i)
+  post_probs   <- apply(contr$contMat, 2, getPostProb, post_combs_i)
   
-  res <- c(sign       = ifelse(max(post_probs) > crit_prob, 1, 0),
-           p_val      = max(post_probs),
-           post_probs = post_probs)
+  res <- c(sign          = ifelse(max(post_probs) > crit_prob_adj, 1, 0),
+           crit_prob_adj = crit_prob_adj,
+           max_post_prob = max(post_probs),
+           post_probs    = post_probs)
   
   return (res)
   
-}
+}
\ No newline at end of file
diff --git a/R/bootstrapping.R b/R/bootstrapping.R
index 8585b83..7570a8e 100644
--- a/R/bootstrapping.R
+++ b/R/bootstrapping.R
@@ -1,20 +1,40 @@
-#' @title getBootstrapBands
+#' @title getBootstrapQuantiles
 #'
-#' @param model_fits tbd
-#' @param n_samples tbd
-#' @param alpha tbd
-#' @param avg_fit tbd
-#' @param dose_seq tbd
+#' @description A function for the calculation of credible intervals to assess the uncertainty for the model fit. 
+#' Hereby the credible intervals are calculated as follows.
+#' Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed. 
+#' These fits are then used to identify the specified quantiles. 
+#' This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley et al. (2017).
+#' Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution.
+#' @references O'Quigley, J., Iasonos, A., & Bornkamp, B. (Eds.). (2017). Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315151984
+#' @param model_fits an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting 
+#' @param quantiles a vector of quantiles that should be evaluated 
+#' @param n_samples number of samples that should be drawn as basis for the 
+#' @param doses a vector of doses for which a prediction should be performed
+#' @param avg_fit boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.
 #'
-#' @return tbd
+#' @examples
+#' # example code
+#' posterior_list =   list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+#'                    DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+#'                    DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+#'                    DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+#'                    DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2))
+#' models=c("emax","exponential","sigEmax","linear")
+#' dose_levels=c(0,1,2,4,8)
+#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels)
+#' fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE)
+#' getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8))
+#' getBootstrapQuantiles(fit_simple, n_samples=2000, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8))
+#' @return  A data frame with entries doses, models, and quantiles 
 #' @export
-getBootstrapBands <- function (
+getBootstrapQuantiles <- function (
 
   model_fits,
+  quantiles,
   n_samples = 1e3,
-  alpha     = c(0.05, 0.5),
-  avg_fit   = TRUE,
-  dose_seq  = NULL
+  doses     = NULL,
+  avg_fit   = TRUE
 
 ) {
   
@@ -24,11 +44,11 @@ getBootstrapBands <- function (
   
   dose_levels    <- model_fits[[1]]$dose_levels
   model_names    <- names(model_fits)
-  quantile_probs <- c(0.5, sort(unique(c(alpha / 2, 1 - alpha / 2))))
+  quantile_probs <- sort(unique(quantiles))
   
-  if (is.null(dose_seq)) {
+  if (is.null(doses)) {
     
-    dose_seq <- seq(min(dose_levels), max(dose_levels), length.out = 100L)
+    doses <- seq(min(dose_levels), max(dose_levels), length.out = 100L)
     
   }
   
@@ -45,7 +65,7 @@ getBootstrapBands <- function (
         bnds  = DoseFinding::defBnds(
           mD = max(model_fits[[1]]$dose_levels))[[model]])
       
-      preds <- stats::predict(fit, doseSeq = dose_seq, predType = "ls-means")
+      preds <- stats::predict(fit, doseSeq = doses, predType = "ls-means")
       attr(preds, "gAIC") <- DoseFinding::gAIC(fit)
       
       return (preds)
@@ -71,20 +91,20 @@ getBootstrapBands <- function (
     
   }
   
-  sort_indx <- order(rep(seq_along(model_names), length(dose_seq)))
+  sort_indx <- order(rep(seq_along(model_names), length(doses)))
   quant_mat <- t(apply(X      = preds[sort_indx, ],
                        MARGIN = 1,
                        FUN    = stats::quantile,
                        probs  = quantile_probs))
   
   cr_bounds_data <- cbind(
-    dose_seqs = dose_seq,
-    models    = rep(
+    doses  = doses,
+    models = rep(
       x    = factor(model_names,
                     levels = c("linear", "emax", "exponential",
                                "sigEmax", "logistic", "quadratic",
                                "avgFit")),
-      each = length(dose_seq)),
+      each = length(doses)),
     as.data.frame(quant_mat))
 
   return (cr_bounds_data)
diff --git a/R/modelling.R b/R/modelling.R
index 1328d24..ed50b13 100644
--- a/R/modelling.R
+++ b/R/modelling.R
@@ -1,15 +1,39 @@
 #' @title getModelFits
 #' 
-#' @description Fits dose-response curves for the specified dose-repsonse models, based on the posterior distributions.
+#' @description Fits dose-response curves for the specified dose-response models, based on the posterior distributions.
 #' For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions. 
-#' For the default case, the Nelder-Mead algorithm is used. Will be further updated and links to publication as well as references will be added.
-#' 
+#' For the default case, the Nelder-Mead algorithm is used. 
+#' In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for \eqn{l \in {1,...,L}} are used as basis
+#' For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m}
+#' \deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{l=1}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))}
+#' Herefore the function nloptr of the nloptr package is utilized. 
+#' In the simplified case \eqn{L=1}, as the dimension of the posterior is reduced to 1 first.
+#' The generalized AIC values are calculated via the formula
+#' \deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p }
+#' where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels.
+#' Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}.
+#' Corresponding gAIC based weights for model \eqn{M} are calculated as outlined in Schorning et al. (2016)
+#' \deqn{
+#' \Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_{m})}
+#' }
+#' where \eqn{Q} denotes the number of models included in the averaging procedure.
+#' @references Schorning K, Bornkamp B, Bretz F, Dette H. 2016. “Model selection versus model averaging in dose finding studies”. Stat Med; 35; 4021-4040.
 #' @param models list of model names for which a fit will be performed.
 #' @param dose_levels a vector containing the different dosage levels.
 #' @param posterior a getPosterior object, containing the (multivariate) posterior distribution per dosage level. 
 #' @param simple boolean variable, defining whether simplified fit will be applied. Default FALSE.
-#' 
-#' @return model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC per model.
+#' @examples
+#' # example code
+#' posterior_list =   list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+#'                    DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+#'                    DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+#'                    DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+#'                    DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2))
+#' models=c("emax","exponential","sigEmax","linear")
+#' dose_levels=c(0,1,2,4,8)
+#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels)
+#' fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE)
+#' @return model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC (and corresponding weight) per model.
 #' 
 #' @export
 getModelFits <- function (
@@ -21,6 +45,8 @@ getModelFits <- function (
   
 ) {
   
+  models      <- unique(gsub("\\d", "", models))
+  
   getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt)
   model_fits  <- lapply(models, getModelFit, dose_levels, posterior)
   
@@ -30,6 +56,7 @@ getModelFits <- function (
   attr(model_fits, "posterior") <- posterior
   class(model_fits)             <- "modelFits"
   
+ 
   return (model_fits)
   
 }
diff --git a/R/plot.R b/R/plot.R
index efd1a4a..b53821c 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -1,7 +1,9 @@
 #' @title plot.modelFits
 #' 
 #' @description plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit.
-#' More details to be added, as well as references.
+#' Black lines show the fitted dose response models and an AIC based average model. Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group). 
+#' To assess the uncertainty of the model fit one can in addition visualize credible bands (default coloring as orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function.
+#' The default setting is that these credible bands are not calculated. 
 #' @param x an object of type getModelFits
 #' @param gAIC logical value indicating whether gAIC values are shown in the plot. Default TRUE
 #' @param avg_fit logical value indicating whether average fit is presented in the plot. Default TRUE
@@ -12,7 +14,19 @@
 #' @param n_bs_smpl number of bootstrap samples being used. Default set to 1000.
 #' @param acc_color color of the credible bands. Default set to "orange"
 #' @param ... optional parameter to be passed.
-#' 
+#' @examples
+#' # example code
+#' posterior_list =   list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+#'                    DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+#'                    DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+#'                    DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+#'                    DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2))
+#' models=c("emax","exponential","sigEmax","linear")
+#' dose_levels=c(0,1,2,4,8)
+#' fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels)
+#' fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE)
+#' plot(fit, cr_bands = TRUE)
+#' plot(fit_simple, cr_bands = TRUE,alpha_CrB =c(0.05,0.1,0.5))
 #' @return plts returns a ggplot2 object
 #' @export
 plot.modelFits <- function (
@@ -37,27 +51,27 @@ plot.modelFits <- function (
   post_summary <- summary.postList(
     object = attr(model_fits, "posterior"),
     probs  = c(alpha_CrI / 2, 0.5, 1 - alpha_CrI / 2))
-  doses        <- seq(from       = min(dose_levels),
-                      to         = max(dose_levels),
-                      length.out = plot_res)
+  dose_seq <- seq(from       = min(dose_levels),
+                  to         = max(dose_levels),
+                  length.out = plot_res)
   
-  preds_models <- sapply(model_fits, predictModelFit, doses = doses)
+  preds_models <- sapply(model_fits, predictModelFit, doses = dose_seq)
   model_names  <- names(model_fits)
   
   if (avg_fit) {
     
     mod_weights <- sapply(model_fits, function (x) x$model_weight)
-    avg_mod     <- preds_models %*% mod_weights
+    avg_preds   <- preds_models %*% mod_weights
     
-    preds_models <- cbind(preds_models, avg_mod)
+    preds_models <- cbind(preds_models, avg_preds)
     model_names  <- c(model_names, "avgFit")
     
   }
   
   gg_data <- data.frame(
-    dose_seqs = rep(doses, length(model_names)),
-    fits      = as.vector(preds_models),
-    models    = rep(factor(model_names,
+    doses  = rep(dose_seq, length(model_names)),
+    fits   = as.vector(preds_models),
+    models = rep(factor(model_names,
                            levels = c("linear", "emax", "exponential",
                                       "sigEmax", "logistic", "quadratic",
                                       "avgFit")),
@@ -88,12 +102,12 @@ plot.modelFits <- function (
   
   if (cr_bands) {
     
-    crB_data <- getBootstrapBands(
+    crB_data <- getBootstrapQuantiles(
       model_fits = model_fits,
       n_samples  = n_bs_smpl,
-      alpha      = alpha_CrB,
+      quantiles  = c(0.5, sort(unique(c(alpha_CrB / 2, 1 - alpha_CrB / 2)))),
       avg_fit    = avg_fit,
-      dose_seq   = doses)
+      doses      = dose_seq)
     
     getInx <- function (alpha_CrB) {
       n        <- length(alpha_CrB)
@@ -131,7 +145,7 @@ plot.modelFits <- function (
       loop_txt <- paste0(
         "ggplot2::geom_ribbon(
           data    = crB_data,
-          mapping = ggplot2::aes(x    = dose_seqs,
+          mapping = ggplot2::aes(x    = doses,
                                  ymin = crB_data[, ", inx[1], "],
                                  ymax = crB_data[, ", inx[2], "]),
           fill    = acc_color,                    
@@ -146,7 +160,7 @@ plot.modelFits <- function (
     plts <- plts +
       ggplot2::geom_line(
         data    = crB_data,
-        mapping = ggplot2::aes(dose_seqs, `50%`),
+        mapping = ggplot2::aes(doses, `50%`),
         color   = acc_color)
     
   }
@@ -177,7 +191,7 @@ plot.modelFits <- function (
     ## Fitted Models
     ggplot2::geom_line(
       data    = gg_data,
-      mapping = ggplot2::aes(dose_seqs, fits)) + 
+      mapping = ggplot2::aes(doses, fits)) + 
     ## Faceting
     ggplot2::facet_wrap(~ models)
   
diff --git a/R/posterior.R b/R/posterior.R
index 5c12d44..6693fbb 100644
--- a/R/posterior.R
+++ b/R/posterior.R
@@ -1,104 +1,61 @@
-#' @title getPriorList
-#' 
-#' @param hist_data historical trial summary level data,
-#' needs to be provided as a dataframe. Including information of the
-#' estimates and variability.
-#' @param dose_levels vector of the different doseage levels
-#' @param dose_names character vector of dose levels,
-#' default NULL and will be automatically created
-#' based on the dose levels parameter.
-#' @param robustify_weight Null needs to be provided as a numeric
-#' value for the weight of the robustification component
-#'
-#' @export
-getPriorList <- function (
-    
-  hist_data,
-  dose_levels,
-  dose_names       = NULL,
-  robustify_weight = NULL
-  
-) {
-  
-  sd_tot <- with(hist_data, sum(sd * n) / sum(n))
-  
-  gmap <- RBesT::gMAP(
-    formula    = cbind(est, se) ~ 1 | trial,
-    weights    = hist_data$n,
-    data       = hist_data,
-    family     = gaussian,
-    beta.prior = cbind(0, 100 * sd_tot),
-    tau.dist   = "HalfNormal",
-    tau.prior  = cbind(0, sd_tot / 4))
-  
-  prior_ctr <- RBesT::automixfit(gmap)
-  
-  if(is.null(robustify_weight) | !is.numeric(robustify_weight)) {
-    stop("robustify_weight needs to be provided and must be numeric")
-  }
-  
-  if (!is.null(robustify_weight)) {
-    
-    prior_ctr <- suppressMessages(RBesT::robustify(
-      priormix = prior_ctr,
-      weight   = robustify_weight,
-      sigma    = sd_tot))
-    
-  }
-  
-  prior_trt <- RBesT::mixnorm(
-    comp1 = c(w = 1, m = summary(prior_ctr)[1], n = 1),
-    sigma = sd_tot,
-    param = "mn")
-  
-  prior_list <- c(list(prior_ctr),
-                  rep(x     = list(prior_trt),
-                      times = length(dose_levels[-1])))
-  
-  if (is.null(dose_names)) {
-    
-    dose_names <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1])))
-    
-  }
-  
-  names(prior_list) <- dose_names
-  
-  attr(prior_list, "dose_levels") <- dose_levels
-  attr(prior_list, "sd_tot")      <- sd_tot
-  
-  return (prior_list)
-  
-}
 
 #' @title getPosterior
 #' 
-#' @description Either the patient level data or both the mu_hat as well as the sd_hat must to be provided.
+#' @description Either the patient level data or both mu_hat as well as sd_hat must to be provided. If patient level data is provided mu_hat and se_hat are calculated within the function using a linear model.
+#' This function calculates the posterior for every dose group independently via the RBesT function postmix.
 #' 
-#' @param data dataframe containing the information of dose and response.
-#' Also a simulateData object can be provided.
 #' @param prior_list prior_list object
-#' @param mu_hat vector of estimated mean values
-#' @param sd_hat vector of estimated standard deviations.
-#'
+#' @param data dataframe containing the information of dose and response. Default NULL
+#' Also a simulateData object can be provided.
+#' @param mu_hat vector of estimated mean values (per dose group).
+#' @param se_hat vector of estimated standard deviations (per dose group).
+#' @param calc_ess boolean variable, indicating whether effective sample size should be calculated. Default FALSE
+#' @return posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group
+#' @examples
+#' # example code
+#' prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
+#'                    DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+#'                    DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+#'                    DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+#'                    DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+#' mu<-c(0,1,1.5,2,2.5)
+#' se<-c(5,4,6,7,8)
+#' posterior_list <- getPosterior(
+#'    prior_list = prior_list,
+#'     mu_hat   = mu,
+#'    se_hat   = se)
+#'    
+#' summary(posterior_list)
+#' 
 #' @export
 getPosterior <- function(
-  data = NULL,
   prior_list,
-  mu_hat = NULL,
-  sd_hat = NULL
+  data     = NULL,
+  mu_hat   = NULL,
+  se_hat   = NULL,
+  calc_ess = FALSE
   
 ) {
-  if (is.null(data)){
-    posterior_list <-getPosteriorI(data_i=NULL, prior_list = prior_list,
-                                   mu_hat     = mu_hat,
-                                   sd_hat     = sd_hat)
-  }else{
-  posterior_list <- lapply(split(data, data$simulation), getPosteriorI,
-                           prior_list = prior_list,
-                           mu_hat     = mu_hat,
-                           sd_hat     = sd_hat)
-  }
   
+  if (!is.null(mu_hat) && !is.null(se_hat) && is.null(data)) {
+    
+    posterior_list <- getPosteriorI(
+      prior_list = prior_list,
+      mu_hat     = mu_hat,
+      se_hat     = se_hat,
+      calc_ess   = calc_ess)
+    
+  } else if (is.null(mu_hat) && is.null(se_hat) && !is.null(data)) {
+    
+    posterior_list <- lapply(split(data, data$simulation), getPosteriorI,
+                             prior_list = prior_list, calc_ess = calc_ess)
+    
+  } else {
+    
+    stop ("Either 'data' or 'mu_hat' and 'se_hat' must not be NULL.")
+    
+  }
+ 
   if (length(posterior_list) == 1) {
     
     posterior_list <- posterior_list[[1]]
@@ -111,33 +68,35 @@ getPosterior <- function(
 
 getPosteriorI <- function(
     
-  data_i,
+  data_i   = NULL,
   prior_list,
-  mu_hat = NULL,
-  sd_hat = NULL
+  mu_hat   = NULL,
+  se_hat   = NULL,
+  calc_ess = FALSE
   
 ) {
   
-  if (is.null(mu_hat) && is.null(sd_hat)) {
+  if (is.null(mu_hat) && is.null(se_hat)) {
     
     anova_res <- stats::lm(data_i$response ~ factor(data_i$dose) - 1)
     mu_hat    <- summary(anova_res)$coefficients[, 1]
-    sd_hat    <- summary(anova_res)$coefficients[, 2]
+    se_hat    <- summary(anova_res)$coefficients[, 2]
     
-  } else if (!is.null(mu_hat) && !is.null(sd_hat)) {
+  } else if (!is.null(mu_hat) && !is.null(se_hat)) {
     
     stopifnot("m_hat length must match number of dose levels" = 
                 length(prior_list) == length(mu_hat),
-              "sd_hat length must match number of dose levels" = 
-                length(prior_list) == length(sd_hat))
+              "se_hat length must match number of dose levels" = 
+                length(prior_list) == length(se_hat))
     
   } else {
     
-    stop ("Both mu_hat and sd_hat must be provided.")
+    stop ("Both mu_hat and se_hat must be provided.")
     
   }
   
-  post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = sd_hat)
+  post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = se_hat,
+                      SIMPLIFY = FALSE)
   
   if (is.null(names(prior_list))) {
     
@@ -148,10 +107,38 @@ getPosteriorI <- function(
   names(post_list) <- names(prior_list)
   class(post_list) <- "postList"
   
+  if (calc_ess) {
+    
+    attr(post_list, "ess") <- getESS(post_list)
+    
+  } else {
+    
+    attr(post_list, "ess") <- numeric(0)
+    
+  }
+  
   return (post_list)
   
 }
 
+#' @title getESS
+#' 
+#' @description This function calculates the effective sample size for every dose group via the RBesT function ess.
+#' 
+#' @param post_list a posterior list object, for which the effective sample size (per dose group) should be calculated
+#'
+#' @return a vector of the effective sample sizes (per dose group) 
+#' @export
+getESS <- function (
+    
+  post_list
+  
+) {
+  
+  suppressMessages(sapply(post_list, RBesT::ess))
+  
+}
+
 getPostCombsI <- function (
     
   posterior_i
diff --git a/R/s3methods.R b/R/s3methods.R
index 14aa663..a507df9 100644
--- a/R/s3methods.R
+++ b/R/s3methods.R
@@ -8,9 +8,11 @@ print.BayesianMCPMod <- function (
   
 ) {
   
-  n_models      <- ncol(x$BayesianMCP) - 2L
-  model_names   <- colnames(x$BayesianMCP)[-c(1, 2)] |>
-    sub(pattern = "post_probs.", replacement = "", x = _)
+  model_names <- colnames(x$BayesianMCP)[
+    grepl("post_probs.", colnames(x$BayesianMCP))] |>
+    sub(pattern = "post_probs.", replacement = "", x = _) |>
+    gsub(pattern = "\\d", replacement = "", x = _) |>
+    unique(x = _)
   
   model_success <- colMeans(do.call(rbind, lapply(x$Mod, function (y) {
     
@@ -20,7 +22,7 @@ print.BayesianMCPMod <- function (
       
     } else {
       
-      model_signs <- rep(FALSE, n_models)
+      model_signs        <- rep(FALSE, length(model_names))
       names(model_signs) <- model_names
       
       return (model_signs)
@@ -34,6 +36,13 @@ print.BayesianMCPMod <- function (
   cat("Model Significance Frequencies\n")
   print(model_success, ...)
   
+  if (any(!is.na(attr(x$BayesianMCP, "ess_avg")))) {
+    
+    cat("Average Posterior ESS\n")
+    print(attr(x$BayesianMCP, "ess_avg"), ...)
+    
+  }
+  
 }
 
 ## BayesianMCP --------------------------------------------
@@ -46,27 +55,65 @@ print.BayesianMCP <- function (
   
 ) {
   
-  power <- mean(x[, 1])
   n_sim <- nrow(x)
   
   cat("Bayesian Multiple Comparison Procedure\n")
-  cat("  Estimated Success Rate: ", power, "\n")
-  cat("  N Simulations:          ", n_sim)
+  
+  if (n_sim == 1L) {
+    
+    attr(x, "crit_prob_adj") <- NULL
+    attr(x, "success_rate")  <- NULL
+    class(x) <- NULL
+    
+    print.default(x, ...)
+    
+  } else {
+    
+    cat("  Estimated Success Rate: ", attr(x, "successRate"), "\n")
+    cat("  N Simulations:          ", n_sim)
+    
+  }
   
 }
 
 ## ModelFits ----------------------------------------------
-
+#' @title predict.modelFits
+#' @description This function performs model predictions based on the provided
+#' model and dose specifications 
+#' 
+#' @param object a modelFits object containing information about the fitted
+#' model coefficients
+#' @param doses a vector specifying the doses for which a prediction should be
+#' done getContrMat object, contrast matrix to be used for the testing step.
+#' @param ... other parameters
+#' @examples
+#' # example code
+#' posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+#'                        DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+#'                        DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+#'                        DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+#'                        DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1), sigma = 2))
+#' models         <- c("emax", "exponential", "sigEmax", "linear")
+#' dose_levels    <- c(0, 1, 2, 4, 8)
+#' fit            <- getModelFits(models      = models,
+#'                                posterior   = posterior_list,
+#'                                dose_levels = dose_levels)
+#' predict(fit, doses = c(0, 1, 3, 4, 6, 8))
+#' 
+#' @return a list with the model predictions for the specified models and doses  
 #' @export
-predict.ModelFits <- function (
+predict.modelFits <- function (
     
   object,
   doses = NULL,
   ...
   
 ) {
+
+  predictions <- lapply(object, predictModelFit, doses = doses)
+  attr(predictions, "doses") <- doses
   
-  lapply(object, predictModelFit, doses = doses, ...)
+  return (predictions)
   
 }
 
@@ -88,7 +135,8 @@ print.modelFits <- function (
   
   out_table <- data.frame(predictions,
                           mEff = sapply(x, function (y) y$max_effect),
-                          gAIC = sapply(x, function (y) y$gAIC))
+                          gAIC = sapply(x, function (y) y$gAIC),
+                          w    = sapply(x, function (y) y$model_weight))
   out_table <- apply(as.matrix(out_table), 2, round, digits = n_digits)
   
   if (!is.null(x[[1]]$significant)) {
@@ -122,7 +170,7 @@ print.modelFits <- function (
   cat("Dose Levels\n",
       paste(dose_names, round(dose_levels, n_digits), sep = " = "), "\n")
   cat("\n")
-  cat("Predictions, Maximum Effect, gAIC & Significance\n")
+  cat("Predictions, Maximum Effect, gAIC, Model Weights & Significance\n")
   print(out_table, ...)
   
 }
diff --git a/R/simulation.R b/R/simulation.R
index 3eb16ab..c4bd480 100644
--- a/R/simulation.R
+++ b/R/simulation.R
@@ -1,16 +1,37 @@
 #' @title simulateData
 #' 
+#' @description
+#' Function to simulate patient level data for a normally distributed endpoint
+#' 
+#' 
 #' @param n_patients vector containing number of patients as a numerical
 #' value per dose-group.
-#' @param dose_levels vector containing the different doseage levels.
+#' @param dose_levels vector containing the different dosage levels.
 #' @param sd standard deviation on patient level.
-#' @param mods An object of class "Mods" as specified in the Dosefinding package.
+#' @param mods An object of class "Mods" as specified in the Doseinding package.
 #' @param n_sim number of simulations to be performed,
 #' Default is 1000
 #' @param true_model Default value is NULL.
 #' Assumed true underlying model. Provided via a String. e.g. "emax".
 #' In case of NULL, all dose-response models, included in the mods input parameter will be used.
+#' @param dr_means a vector, with information about assumed effects per dose group. Default NULL.
 #' 
+#' @examples
+#' # example code
+#'  models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, 
+#' doses = c(0, 0.5, 2,4, 8),maxEff= 6)
+#' dose_levels = c(0, 0.5, 2,4, 8)
+#' sd   = 12
+#' n_patients <- c(40,60,60,60,60)
+#' sim_data <- simulateData(
+#' n_patients  = n_patients,
+#' dose_levels = dose_levels,
+#' mods        = models,
+#' n_sim       = 100,
+#' sd          = sd)
+#' 
+#' sim_data
+
 #' @return sim_data one list object, containing patient level simulated data for all assumed true models.
 #' Also providing information about simulation iteration, patient number as well as dosage levels.
 #' 
@@ -21,8 +42,9 @@ simulateData <- function(
   dose_levels,
   sd,
   mods,
-  n_sim = 1e3,
-  true_model = NULL
+  n_sim      = 1e3,
+  true_model = NULL,
+  dr_means   = NULL
   
 ) {
   
@@ -41,10 +63,24 @@ simulateData <- function(
     ptno       = rep(seq_len(sum(n_patients)), times = n_sim),
     dose       = rep(rep(dose_levels, times = n_patients), times = n_sim))
   
-  model_responses <- DoseFinding::getResp(mods, sim_info$dose)
-  random_noise    <- stats::rnorm(nrow(sim_info), mean = 0, sd = sd)
+  if (is.null(dr_means)) {
+    
+    model_responses <- DoseFinding::getResp(mods, sim_info$dose)
+    
+  } else {
+    
+    stopifnot(identical(length(dose_levels), length(dr_means)))
+    
+    model_responses <- matrix(
+      data     = rep(unlist(mapply(rep, dr_means, n_patients)),
+                     n_sim * length(mods)),
+      ncol     = length(mods),
+      dimnames = list(NULL, names(mods)))
+    
+  }
   
-  sim_data <- cbind(sim_info, model_responses + random_noise)
+  random_noise <- stats::rnorm(nrow(sim_info), mean = 0, sd = sd)
+  sim_data     <- cbind(sim_info, model_responses + random_noise)
   
   if (!is.null(true_model)) {
     
@@ -52,6 +88,13 @@ simulateData <- function(
     
   }
   
+  if (!is.null(dr_means)) {
+    
+    sim_data           <- sim_data[1:4]
+    colnames(sim_data) <- c(colnames(sim_data)[-4], "dose_resp")
+    
+  }
+  
   return (sim_data)
   
 }
@@ -68,4 +111,4 @@ getModelData <- function (
   
   return (model_data)
   
-}
\ No newline at end of file
+}
diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd
index 6883d23..1f0a539 100644
--- a/man/assessDesign.Rd
+++ b/man/assessDesign.Rd
@@ -2,30 +2,65 @@
 % Please edit documentation in R/BMCPMod.R
 \name{assessDesign}
 \alias{assessDesign}
-\title{assessDesign}
+\title{assessDesign
+.}
 \usage{
 assessDesign(
   n_patients,
   mods,
   prior_list,
+  sd,
   n_sim = 1000,
   alpha_crit_val = 0.05,
-  simple = TRUE
+  simple = TRUE,
+  reestimate = FALSE,
+  contr = NULL,
+  dr_means = NULL
 )
 }
 \arguments{
-\item{n_patients}{tbd}
+\item{n_patients}{Vector specifying the planned number of patients per dose group}
 
-\item{mods}{tbd}
+\item{mods}{An object of class "Mods" as specified in the DoseFinding package.}
 
-\item{prior_list}{tbd}
+\item{prior_list}{a prior_list object specifying the utilized prior for the different dose groups}
 
-\item{n_sim}{tbd}
+\item{sd}{a positive value, specification of assumed sd}
 
-\item{alpha_crit_val}{tbd}
+\item{n_sim}{number of simulations to be performed}
 
-\item{simple}{tbd}
+\item{alpha_crit_val}{(unadjusted) critical value to be used for the MCT testing step. Passed to the getCritProb function for the calculation of adjusted critical values (on the probability scale). Default is 0.05.}
+
+\item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.}
+
+\item{reestimate}{boolean variable, defining whether critical value should be calculated with re-estimated contrasts (see getCritProb function for more details). Default FALSE}
+
+\item{contr}{Allows specification of a fixed contrasts matrix. Default NULL}
+
+\item{dr_means}{a vector, allows specification of  individual (not model based) assumed effects per dose group. Default NULL}
+}
+\value{
+returns success probabilities for the different assumed dose-response shapes, attributes also includes information around average success rate (across all assumed models) and prior Effective sample size
 }
 \description{
-assessDesign
+This function performs simulation based trial design evaluations for a set of specified dose-response models
+}
+\examples{
+# example code
+mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6)
+sd   = 12
+prior_list<-list(Ctrl   = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma = 2),
+                   DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+                   DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+                   DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+                   DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+
+n_patients <- c(40,60,60,60,60)
+success_probabilities <- assessDesign(
+n_patients  = n_patients,
+mods        = mods,
+prior_list  = prior_list,
+sd          = sd)
+
+success_probabilities
 }
diff --git a/man/getBootstrapBands.Rd b/man/getBootstrapBands.Rd
deleted file mode 100644
index de4bd80..0000000
--- a/man/getBootstrapBands.Rd
+++ /dev/null
@@ -1,31 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/bootstrapping.R
-\name{getBootstrapBands}
-\alias{getBootstrapBands}
-\title{getBootstrapBands}
-\usage{
-getBootstrapBands(
-  model_fits,
-  n_samples = 1000,
-  alpha = c(0.05, 0.5),
-  avg_fit = TRUE,
-  dose_seq = NULL
-)
-}
-\arguments{
-\item{model_fits}{tbd}
-
-\item{n_samples}{tbd}
-
-\item{alpha}{tbd}
-
-\item{avg_fit}{tbd}
-
-\item{dose_seq}{tbd}
-}
-\value{
-tbd
-}
-\description{
-getBootstrapBands
-}
diff --git a/man/getBootstrapQuantiles.Rd b/man/getBootstrapQuantiles.Rd
new file mode 100644
index 0000000..4e63309
--- /dev/null
+++ b/man/getBootstrapQuantiles.Rd
@@ -0,0 +1,53 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/bootstrapping.R
+\name{getBootstrapQuantiles}
+\alias{getBootstrapQuantiles}
+\title{getBootstrapQuantiles}
+\usage{
+getBootstrapQuantiles(
+  model_fits,
+  quantiles,
+  n_samples = 1000,
+  doses = NULL,
+  avg_fit = TRUE
+)
+}
+\arguments{
+\item{model_fits}{an object of class modelFits, i.e. information about fitted models & corresponding model coefficients as well as the posterior distribution that was the basis for the model fitting}
+
+\item{quantiles}{a vector of quantiles that should be evaluated}
+
+\item{n_samples}{number of samples that should be drawn as basis for the}
+
+\item{doses}{a vector of doses for which a prediction should be performed}
+
+\item{avg_fit}{boolean variable, defining whether an average fit (based on generalized AIC weights) should be performed in addition to the individual models. Default TRUE.}
+}
+\value{
+A data frame with entries doses, models, and quantiles
+}
+\description{
+A function for the calculation of credible intervals to assess the uncertainty for the model fit.
+Hereby the credible intervals are calculated as follows.
+Samples from the posterior distribution are drawn (via the RBesT function rmix) and for every sample the simplified fitting step (see getModelFits function) and a prediction is performed.
+These fits are then used to identify the specified quantiles.
+This approach can be considered as the bayesian equivalent of the frequentist bootstrap approach described in O'Quigley et al. (2017).
+Instead of drawing n bootstrap samples from the sampling distribution of the trial dose-response estimates, here the samples are directly taken from the posterior distribution.
+}
+\examples{
+# example code
+posterior_list =   list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+                   DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+                   DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+                   DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+                   DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2))
+models=c("emax","exponential","sigEmax","linear")
+dose_levels=c(0,1,2,4,8)
+fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels)
+fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE)
+getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8))
+getBootstrapQuantiles(fit_simple, n_samples=2000, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8))
+}
+\references{
+O'Quigley, J., Iasonos, A., & Bornkamp, B. (Eds.). (2017). Handbook of Methods for Designing, Monitoring, and Analyzing Dose-Finding Trials (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315151984
+}
diff --git a/man/getContr.Rd b/man/getContr.Rd
new file mode 100644
index 0000000..c1e5f20
--- /dev/null
+++ b/man/getContr.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/BMCPMod.R
+\name{getContr}
+\alias{getContr}
+\title{getContr}
+\usage{
+getContr(
+  mods,
+  dose_levels,
+  dose_weights = NULL,
+  prior_list = NULL,
+  sd_posterior = NULL,
+  se_new_trial = NULL
+)
+}
+\arguments{
+\item{mods}{An object of class "Mods" as specified in the DoseFinding package.}
+
+\item{dose_levels}{vector containing the different dosage levels.}
+
+\item{dose_weights}{Vector specifying weights for the different doses. Please note that in case this information is provided together with a prior (i.e. Option i) is planned) these two inputs should be provided on the same scale (e.g. patient numbers).  Default NULL}
+
+\item{prior_list}{a prior_list object, only required as input for Option i). Default NULL}
+
+\item{sd_posterior}{a vector of positive values with information about the variability of the posterior distribution, only required for Option iii). Default NULL}
+
+\item{se_new_trial}{a vector of positive values with information about the observed variability, only required for Option iv). Default NULL}
+}
+\value{
+contr Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the DoseFinding package.
+}
+\description{
+This function calculates contrast vectors that are optimal for detecting certain alternatives via applying the function optContr of the DoseFinding package.
+Hereby 4 different options can be distinguished that are automatically executed based on the input that is provided
+i)  Bayesian approach: If dose_weights and a prior_list are provided an optimized contrasts for the posterior sample size is calculated.
+In detail,  in a first step the dose_weights (typically the number of patients per dose group) and the prior information is combined by calculating for
+each dose group a posterior effective sample. Based on this posterior effective sample sizes the allocation ratio is derived, which allows for a calculation on
+pseudo-optimal contrasts via regular MCPMod.are calculated from the
+regular MCPMod for these specific weights
+ii) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
+regular MCPMod for these specific weights
+iii)Bayesian approach + re-estimation: If only a sd_posterior (i.e. variability of the posterior distribution) is provided, pseudo-optimal contrasts based on these posterior weights will be calculated
+iv) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the
+regular MCPMod for this specific vector of standard errors. For the actual evaluation this vector of standard errors is translated into a (diagonal) matrix of variances
+}
+\examples{
+# example code
+mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8), maxEff= 6)
+dose_levels=c(0, 0.5, 2, 4, 8)
+sd_posterior   = c(2.8,3,2.5,3.5,4)
+contr_mat<- getContr(
+mods           = mods,
+dose_levels    = dose_levels,
+sd_posterior   = sd_posterior) 
+
+}
diff --git a/man/getContrMat.Rd b/man/getContrMat.Rd
deleted file mode 100644
index 6040050..0000000
--- a/man/getContrMat.Rd
+++ /dev/null
@@ -1,23 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/BMCPMod.R
-\name{getContrMat}
-\alias{getContrMat}
-\title{getContrMat}
-\usage{
-getContrMat(mods, dose_levels, dose_weights, prior_list)
-}
-\arguments{
-\item{mods}{An object of class "Mods" as specified in the Dosefinding package.}
-
-\item{dose_levels}{vector containing the different doseage levels.}
-
-\item{dose_weights}{Vector specifying weights for the different doses}
-
-\item{prior_list}{a prior_list object}
-}
-\value{
-contr_mat Object of class ‘⁠optContr⁠’. A list containing entries contMat and muMat, and CorrMat. Specified in the Dosefinding package.
-}
-\description{
-This function calculates contrast vectors that are optimal for detecting certain alternatives. More information and link to publication will be added.
-}
diff --git a/man/getCritProb.Rd b/man/getCritProb.Rd
index 21d3482..39593f6 100644
--- a/man/getCritProb.Rd
+++ b/man/getCritProb.Rd
@@ -4,14 +4,22 @@
 \alias{getCritProb}
 \title{getCritProb}
 \usage{
-getCritProb(mods, dose_levels, dose_weights, alpha_crit_val = 0.025)
+getCritProb(
+  mods,
+  dose_levels,
+  dose_weights = NULL,
+  se_new_trial = NULL,
+  alpha_crit_val = 0.025
+)
 }
 \arguments{
-\item{mods}{An object of class "Mods" as specified in the Dosefinding package.}
+\item{mods}{An object of class "Mods" as specified in the DoseFinding package.}
 
-\item{dose_levels}{vector containing the different doseage levels.}
+\item{dose_levels}{vector containing the different dosage levels.}
 
-\item{dose_weights}{Vector specifying weights for the different doses}
+\item{dose_weights}{Vector specifying weights for the different doses, only required for Option i). Default NULL}
+
+\item{se_new_trial}{a vector of positive values, only required for Option ii). Default NULL}
 
 \item{alpha_crit_val}{significance level. Default set to 0.025.}
 }
@@ -19,5 +27,21 @@ getCritProb(mods, dose_levels, dose_weights, alpha_crit_val = 0.025)
 crit_pval multiplicity adjusted critical value on the probability scale.
 }
 \description{
-getCritProb
+This function calculates multiplicity adjusted critical values. The critical values are calculated in such a way that
+when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level).
+Hereby optimal contrasts of the frequentist MCPMod are applied and two options can be distinguished
+i) Frequentist approach: If only dose_weights are provided optimal contrast vectors are calculated from the
+regular MCPMod for these specific weights and the corresponding critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package.
+ii) Frequentist approach+re-estimation:If only a se_new_trial (i.e. the estimated variability per dose group of a new trial) is provided, optimal contrast vectors are calculated from the
+regular MCPMod for this specific vector of standard errors. Here as well the critical value for this set of contrasts is calculated via the critVal function of the DoseFinding package.
+}
+\examples{
+# example code
+mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8))
+dose_levels=c(0, 0.5, 2, 4, 8)
+critVal<- getCritProb(
+  mods           = mods,
+  dose_weights  =c(50,50,50,50,50), #reflecting the planned sample size
+  dose_levels    = dose_levels,
+  alpha_crit_val = 0.05) 
 }
diff --git a/man/getESS.Rd b/man/getESS.Rd
new file mode 100644
index 0000000..7f08776
--- /dev/null
+++ b/man/getESS.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/posterior.R
+\name{getESS}
+\alias{getESS}
+\title{getESS}
+\usage{
+getESS(post_list)
+}
+\arguments{
+\item{post_list}{a posterior list object, for which the effective sample size (per dose group) should be calculated}
+}
+\value{
+a vector of the effective sample sizes (per dose group)
+}
+\description{
+This function calculates the effective sample size for every dose group via the RBesT function ess.
+}
diff --git a/man/getModelFits.Rd b/man/getModelFits.Rd
index 4a87f02..f1b0953 100644
--- a/man/getModelFits.Rd
+++ b/man/getModelFits.Rd
@@ -16,10 +16,39 @@ getModelFits(models, dose_levels, posterior, simple = FALSE)
 \item{simple}{boolean variable, defining whether simplified fit will be applied. Default FALSE.}
 }
 \value{
-model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC per model.
+model_fits returns a list, containing information about the fitted model coefficients, the prediction per dose group as well as maximum effect and generalized AIC (and corresponding weight) per model.
 }
 \description{
-Fits dose-response curves for the specified dose-repsonse models, based on the posterior distributions.
+Fits dose-response curves for the specified dose-response models, based on the posterior distributions.
 For the simplified fit, multivariate normal distributions will be approximated and reduced by one-dimensional normal distributions.
-For the default case, the Nelder-Mead algorithm is used. Will be further updated and links to publication as well as references will be added.
+For the default case, the Nelder-Mead algorithm is used.
+In detail, for both approaches the mean vector \eqn{\theta^{Y}} and the covariance \eqn{\Sigma} of the (mixture) posterior distributions and the corresponding posterior weights \eqn{\tilde{\omega}_{l}} for \eqn{l \in {1,...,L}} are used as basis
+For the full fit a GLS estimator is used to minimize the following expression for the respective dose-response models \eqn{m}
+\deqn{ \hat{\theta}_{m}=argmin_{\theta_{m}} \sum_{l=1}^{L} \tilde{\omega}_{l}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))'\Sigma_{l}^{-1}(\theta_{l_{i}}^{Y}-f(dose_{i},\hat{\theta}_{m}))}
+Herefore the function nloptr of the nloptr package is utilized.
+In the simplified case \eqn{L=1}, as the dimension of the posterior is reduced to 1 first.
+The generalized AIC values are calculated via the formula
+\deqn{gAIC_{m} = \sum_{l=1}^{L} \tilde{\omega}_{l} \sum_{i=0}^{K} \frac{1}{\Sigma_{l_{i,i}}} (\theta_{l_i}^Y - f(dose_{i},\hat{\theta}_{m}))^2 + 2p }
+where \eqn{p} denotes the number of estimated parameters and \eqn{K} the number of active dose levels.
+Here as well for the simplified case the formula reduces to one summand as \eqn{L=1}.
+Corresponding gAIC based weights for model \eqn{M} are calculated as outlined in Schorning et al. (2016)
+\deqn{
+\Omega_I (M) = \frac{\exp(-0.5 gAIC_{M})}{\sum_{m=1}^{Q} \exp(-0.5 gAIC_{m})}
+}
+where \eqn{Q} denotes the number of models included in the averaging procedure.
+}
+\examples{
+# example code
+posterior_list =   list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+                   DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+                   DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+                   DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+                   DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2))
+models=c("emax","exponential","sigEmax","linear")
+dose_levels=c(0,1,2,4,8)
+fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels)
+fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE)
+}
+\references{
+Schorning K, Bornkamp B, Bretz F, Dette H. 2016. “Model selection versus model averaging in dose finding studies”. Stat Med; 35; 4021-4040.
 }
diff --git a/man/getPosterior.Rd b/man/getPosterior.Rd
index 22ac81f..e30d5d2 100644
--- a/man/getPosterior.Rd
+++ b/man/getPosterior.Rd
@@ -4,18 +4,47 @@
 \alias{getPosterior}
 \title{getPosterior}
 \usage{
-getPosterior(data, prior_list, mu_hat = NULL, sd_hat = NULL)
+getPosterior(
+  prior_list,
+  data = NULL,
+  mu_hat = NULL,
+  se_hat = NULL,
+  calc_ess = FALSE
+)
 }
 \arguments{
-\item{data}{dataframe containing the information of dose and response.
+\item{prior_list}{prior_list object}
+
+\item{data}{dataframe containing the information of dose and response. Default NULL
 Also a simulateData object can be provided.}
 
-\item{prior_list}{prior_list object}
+\item{mu_hat}{vector of estimated mean values (per dose group).}
 
-\item{mu_hat}{vector of estimated mean values}
+\item{se_hat}{vector of estimated standard deviations (per dose group).}
 
-\item{sd_hat}{vector of estimated standard deviations.}
+\item{calc_ess}{boolean variable, indicating whether effective sample size should be calculated. Default FALSE}
+}
+\value{
+posterior_list, a posterior list object is returned with information about (mixture) posterior distribution per dose group
 }
 \description{
-Either the patient level data or both the mu_hat as well as the sd_hat must to be provided.
+Either the patient level data or both mu_hat as well as sd_hat must to be provided. If patient level data is provided mu_hat and se_hat are calculated within the function using a linear model.
+This function calculates the posterior for every dose group independently via the RBesT function postmix.
+}
+\examples{
+# example code
+prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
+                   DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+                   DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+                   DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+                   DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+mu<-c(0,1,1.5,2,2.5)
+se<-c(5,4,6,7,8)
+posterior_list <- getPosterior(
+   prior_list = prior_list,
+    mu_hat   = mu,
+   se_hat   = se)
+   
+summary(posterior_list)
+
 }
diff --git a/man/getPriorList.Rd b/man/getPriorList.Rd
deleted file mode 100644
index 7fe0937..0000000
--- a/man/getPriorList.Rd
+++ /dev/null
@@ -1,30 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/posterior.R
-\name{getPriorList}
-\alias{getPriorList}
-\title{getPriorList}
-\usage{
-getPriorList(
-  hist_data,
-  dose_levels,
-  dose_names = NULL,
-  robustify_weight = NULL
-)
-}
-\arguments{
-\item{hist_data}{historical trial summary level data,
-needs to be provided as a dataframe. Including information of the
-estimates and variability.}
-
-\item{dose_levels}{vector of the different doseage levels}
-
-\item{dose_names}{character vector of dose levels,
-default NULL and will be automatically created
-based on the dose levels parameter.}
-
-\item{robustify_weight}{Null needs to be provided as a numeric
-value for the weight of the robustification component}
-}
-\description{
-getPriorList
-}
diff --git a/man/performBayesianMCP.Rd b/man/performBayesianMCP.Rd
index 66e13e5..cc0eaa6 100644
--- a/man/performBayesianMCP.Rd
+++ b/man/performBayesianMCP.Rd
@@ -2,20 +2,53 @@
 % Please edit documentation in R/BMCPMod.R
 \name{performBayesianMCP}
 \alias{performBayesianMCP}
-\title{BayesianMCP}
+\title{performBayesianMCP}
 \usage{
-performBayesianMCP(posteriors_list, contr_mat, crit_prob)
+performBayesianMCP(posterior_list, contr, crit_prob_adj)
 }
 \arguments{
-\item{posteriors_list}{a getPosterior object}
+\item{posterior_list}{a getPosterior object with information about the (mixture) posterior distribution per dose group}
 
-\item{contr_mat}{a getContrMat object, contrast matrix to be used for the testing step.}
+\item{contr}{a getContrMat object, contrast matrix to be used for the testing step.}
 
-\item{crit_prob}{a getCritProb object}
+\item{crit_prob_adj}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale)}
 }
 \value{
-b_mcp test result
+b_mcp test result, with information about p-values for the individual dose-response shapes and overall significance
 }
 \description{
-performs bayesian MCP Test step.
+performs bayesian MCP Test step, as described in Fleischer et al. (2022).
+Tests for a dose-response effect using a model-based multiple contrast test based on the (provided) posterior distribution. In particular for every dose-response candidate the posterior probability is calculated that the contrast is bigger than 0 (based on the posterior distribution of the dose groups).
+In order to obtain significant test decision we consider the maximum of the posterior probabilities across the different models. This maximum is compared with a (multiplicity adjusted) critical value (on the probability scale).
+}
+\examples{
+# example code
+mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8))
+dose_levels=c(0, 0.5, 2, 4, 8)
+sd_posterior   = c(2.8,3,2.5,3.5,4)
+contr_mat<- getContr(
+mods           = mods,
+dose_levels    = dose_levels,
+sd_posterior   = sd_posterior)
+critVal<- getCritProb(
+  mods           = mods,
+  dose_weights  =c(50,50,50,50,50), #reflecting the planned sample size
+  dose_levels    = dose_levels,
+  alpha_crit_val = 0.05)
+prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
+                   DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+                   DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+                   DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+                   DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+mu<-c(0,1,1.5,2,2.5)
+se<-c(5,4,6,7,8)
+posterior_list <- getPosterior(
+   prior_list = prior_list,
+    mu_hat   = mu,
+   se_hat   = se)
+performBayesianMCP(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal)
+
+}
+\references{
+Fleischer F, Bossert S, Deng Q, Loley C, Gierse J. Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670. doi:10.1002/pst.2193
 }
diff --git a/man/performBayesianMCPMod.Rd b/man/performBayesianMCPMod.Rd
index 1bc33dc..b8745ef 100644
--- a/man/performBayesianMCPMod.Rd
+++ b/man/performBayesianMCPMod.Rd
@@ -4,14 +4,14 @@
 \alias{performBayesianMCPMod}
 \title{performBayesianMCPMod}
 \usage{
-performBayesianMCPMod(posteriors_list, contr_mat, crit_prob, simple = FALSE)
+performBayesianMCPMod(posterior_list, contr, crit_prob_adj, simple = FALSE)
 }
 \arguments{
-\item{posteriors_list}{a getPosterior object}
+\item{posterior_list}{a getPosterior object with information about the (mixture) posterior distribution per dose group}
 
-\item{contr_mat}{a getContrMat object, contrast matrix to be used for the testing step.}
+\item{contr}{a getContrMat object, contrast matrix to be used for the testing step.}
 
-\item{crit_prob}{a getCritProb object}
+\item{crit_prob_adj}{a getCritProb object, specifying the critical value to be used for the testing (on the probability scale).}
 
 \item{simple}{boolean variable, defining whether simplified fit will be applied. Passed to the getModelFits function. Default FALSE.}
 }
@@ -19,5 +19,33 @@ performBayesianMCPMod(posteriors_list, contr_mat, crit_prob, simple = FALSE)
 bmcpmod test result as well as modelling result.
 }
 \description{
-performs bayesian MCP Test step and modelling.
+performs bayesian MCP Test step and modelling in a combined fashion. See performBayesianMCP function for MCT Test step and getModelFits for the modelling step
+}
+\examples{
+# example code
+mods <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, doses = c(0, 0.5, 2,4, 8))
+dose_levels=c(0, 0.5, 2, 4, 8)
+sd_posterior   = c(2.8,3,2.5,3.5,4)
+contr_mat<- getContr(
+mods           = mods,
+dose_levels    = dose_levels,
+sd_posterior   = sd_posterior)
+critVal<- getCritProb(
+  mods           = mods,
+  dose_weights  =c(50,50,50,50,50), #reflecting the planned sample size
+  dose_levels    = dose_levels,
+  alpha_crit_val = 0.05)
+prior_list<-list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 5), sigma = 2),
+                   DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 1, s = 12), sigma = 2),
+                   DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 1.2, s = 11), sigma = 2) ,  
+                   DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 1.3, s = 11), sigma = 2) ,
+                   DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 2, s = 13) ,sigma = 2))
+mu<-c(0,1,1.5,2,2.5)
+se<-c(5,4,6,7,8)
+posterior_list <- getPosterior(
+   prior_list = prior_list,
+    mu_hat   = mu,
+   se_hat   = se)
+performBayesianMCPMod(posterior_list=posterior_list, contr=contr_mat,crit_prob_adj=critVal,simple = FALSE)
+
 }
diff --git a/man/plot.modelFits.Rd b/man/plot.modelFits.Rd
index 0cb1c26..8ec1128 100644
--- a/man/plot.modelFits.Rd
+++ b/man/plot.modelFits.Rd
@@ -43,5 +43,21 @@ plts returns a ggplot2 object
 }
 \description{
 plot function based on the ggplot2 package. Providing visualizations for each model and a average Fit.
-More details to be added, as well as references.
+Black lines show the fitted dose response models and an AIC based average model. Dots indicate the posterior median and vertical lines show corresponding credible intervals (i.e. the variability of the posterior distribution of the respective dose group).
+To assess the uncertainty of the model fit one can in addition visualize credible bands (default coloring as orange shaded areas). The calculation of these bands is performed via the getBootstrapQuantiles function.
+The default setting is that these credible bands are not calculated.
+}
+\examples{
+# example code
+posterior_list =   list(Ctrl=RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+                   DG_1=RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+                   DG_2=RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+                   DG_3=RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+                   DG_4=RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1) ,sigma = 2))
+models=c("emax","exponential","sigEmax","linear")
+dose_levels=c(0,1,2,4,8)
+fit<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels)
+fit_simple<-getModelFits(models=models, posterior=posterior_list,dose_levels=dose_levels,simple=TRUE)
+plot(fit, cr_bands = TRUE)
+plot(fit_simple, cr_bands = TRUE,alpha_CrB =c(0.05,0.1,0.5))
 }
diff --git a/man/predict.modelFits.Rd b/man/predict.modelFits.Rd
new file mode 100644
index 0000000..7563063
--- /dev/null
+++ b/man/predict.modelFits.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/s3methods.R
+\name{predict.modelFits}
+\alias{predict.modelFits}
+\title{predict.modelFits}
+\usage{
+\method{predict}{modelFits}(object, doses = NULL, ...)
+}
+\arguments{
+\item{object}{a modelFits object containing information about the fitted
+model coefficients}
+
+\item{doses}{a vector specifying the doses for which a prediction should be
+done getContrMat object, contrast matrix to be used for the testing step.}
+
+\item{...}{other parameters}
+}
+\value{
+a list with the model predictions for the specified models and doses
+}
+\description{
+This function performs model predictions based on the provided
+model and dose specifications
+}
+\examples{
+# example code
+posterior_list <- list(Ctrl = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 1), sigma = 2),
+                       DG_1 = RBesT::mixnorm(comp1 = c(w = 1, m = 3, s = 1.2), sigma = 2),
+                       DG_2 = RBesT::mixnorm(comp1 = c(w = 1, m = 4, s = 1.5), sigma = 2) ,  
+                       DG_3 = RBesT::mixnorm(comp1 = c(w = 1, m = 6, s = 1.2), sigma = 2) ,
+                       DG_4 = RBesT::mixnorm(comp1 = c(w = 1, m = 6.5, s = 1.1), sigma = 2))
+models         <- c("emax", "exponential", "sigEmax", "linear")
+dose_levels    <- c(0, 1, 2, 4, 8)
+fit            <- getModelFits(models      = models,
+                               posterior   = posterior_list,
+                               dose_levels = dose_levels)
+predict(fit, doses = c(0, 1, 3, 4, 6, 8))
+
+}
diff --git a/man/simulateData.Rd b/man/simulateData.Rd
index 9afd8da..86db979 100644
--- a/man/simulateData.Rd
+++ b/man/simulateData.Rd
@@ -10,18 +10,19 @@ simulateData(
   sd,
   mods,
   n_sim = 1000,
-  true_model = NULL
+  true_model = NULL,
+  dr_means = NULL
 )
 }
 \arguments{
 \item{n_patients}{vector containing number of patients as a numerical
 value per dose-group.}
 
-\item{dose_levels}{vector containing the different doseage levels.}
+\item{dose_levels}{vector containing the different dosage levels.}
 
 \item{sd}{standard deviation on patient level.}
 
-\item{mods}{An object of class "Mods" as specified in the Dosefinding package.}
+\item{mods}{An object of class "Mods" as specified in the Doseinding package.}
 
 \item{n_sim}{number of simulations to be performed,
 Default is 1000}
@@ -29,11 +30,29 @@ Default is 1000}
 \item{true_model}{Default value is NULL.
 Assumed true underlying model. Provided via a String. e.g. "emax".
 In case of NULL, all dose-response models, included in the mods input parameter will be used.}
+
+\item{dr_means}{a vector, with information about assumed effects per dose group. Default NULL.}
 }
 \value{
 sim_data one list object, containing patient level simulated data for all assumed true models.
 Also providing information about simulation iteration, patient number as well as dosage levels.
 }
 \description{
-simulateData
+Function to simulate patient level data for a normally distributed endpoint
+}
+\examples{
+# example code
+ models <- DoseFinding::Mods(linear = NULL, linlog = NULL, emax = c(0.5, 1.2), exponential = 2, 
+doses = c(0, 0.5, 2,4, 8),maxEff= 6)
+dose_levels = c(0, 0.5, 2,4, 8)
+sd   = 12
+n_patients <- c(40,60,60,60,60)
+sim_data <- simulateData(
+n_patients  = n_patients,
+dose_levels = dose_levels,
+mods        = models,
+n_sim       = 100,
+sd          = sd)
+
+sim_data
 }
diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R
deleted file mode 100644
index c4caf87..0000000
--- a/tests/testthat/test-bootstrapping.R
+++ /dev/null
@@ -1,152 +0,0 @@
-# Test cases
-test_that("test getBootstrapBands", {
-  library(clinDR)
-  library(dplyr)
-
-  data("metaData")
-  testdata    <- as.data.frame(metaData)
-  dataset     <- filter(testdata, bname == "VIAGRA")
-  histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "ERECTILE DYSFUNCTION")
-
-  ##Create MAP Prior
-  hist_data <- data.frame(
-    trial = c(1, 2, 3, 4),
-    est   = histcontrol$rslt,
-    se    = histcontrol$se,
-    sd    = histcontrol$sd,
-    n     = histcontrol$sampsize)
-
-  sd_tot <- with(hist_data, sum(sd * n) / sum(n))
-
-
-  gmap <- RBesT::gMAP(
-    formula    = cbind(est, se) ~ 1 | trial,
-    weights    = hist_data$n,
-    data       = hist_data,
-    family     = gaussian,
-    beta.prior = cbind(0, 100 * sd_tot),
-    tau.dist   = "HalfNormal",
-    tau.prior  = cbind(0, sd_tot / 4))
-
-  prior_ctr <- RBesT::robustify(
-    priormix = RBesT::automixfit(gmap),
-    weight   = 0.5,
-    mean     = 1.4,
-    sigma    = sd_tot)
-
-  #RBesT::ess(prior_ctr)
-
-  ## derive prior for treatment
-  ## weak informative using same parameters as for robustify component
-  prior_trt <- RBesT::mixnorm(
-    comp1 = c(w = 1, m = 1.4, n = 1),
-    sigma = sd_tot,
-    param = "mn")
-  dose_levels <- c(0, 50, 100, 200)
-  ## combine priors in list
-  prior_list <- c(list(prior_ctr), rep(list(prior_trt), times = length(dose_levels[-1])))
-
-  #Pre-Specification (B)MCPMod
-
-  ## candidate models for MCPMod
-  # linear function - no guestimates needed
-  exp <- DoseFinding::guesst(d     = 50,
-                             p     = c(0.2),
-                             model = "exponential",
-                             Maxd  = max(dose_levels))
-  emax <- DoseFinding::guesst(d     = 100,
-                              p     = c(0.9),
-                              model = "emax")
-
-
-  mods <- DoseFinding::Mods(
-    linear      = NULL,
-    emax        = emax,
-    exponential = exp,
-    doses       = dose_levels,
-    maxEff      = 10,
-    placEff     = 1.4)
-
-  #Simulation of new trial
-  ##Note: This part will be simplified and direct results from one trial will be used
-  mods_sim <- DoseFinding::Mods(
-    emax        = emax,
-    doses       = dose_levels,
-    maxEff      = 12,
-    placEff     = 1.4)
-
-  n_patients <- c(50, 50, 50, 50)
-  data <- simulateData(
-    n_patients  = n_patients,
-    dose_levels = dose_levels,
-    sd          = sd_tot,
-    mods        = mods_sim,
-    n_sim       = 1)
-
-  data_emax           <- data[, c("simulation", "dose", "emax")]
-  names(data_emax)[3] <- "response"
-
-  posterior_emax <- getPosterior(
-    data       = data_emax,
-    prior_list = prior_list)
-
-  #Evaluation of Bayesian MCPMod
-
-  contr_mat <- DoseFinding::optContr(
-    models = mods,
-    doses  = dose_levels,
-    w      = n_patients)
-  ##Calculation of critical value can be done with critVal function
-  crit_val_equal <- DoseFinding:::critVal(contr_mat$corMat, alpha = 0.05, df = 0, alternative = "one.sided")
-  crit_pval      <- pnorm(crit_val_equal)
-
-  ess_prior <- round(unlist(lapply(prior_list, RBesT::ess)))
-
-  ### Evaluation of Bayesian MCPMod
-  contr_mat_prior <- DoseFinding::optContr(
-    models = mods,
-    doses  = dose_levels,
-    w      = n_patients + ess_prior)
-
-  BMCP_result <- performBayesianMCP(
-    posteriors_list = list(posterior_emax),
-    contr_mat       = contr_mat_prior,
-    crit_prob       = crit_pval)
-
-  BMCP_result
-
-  #Model fit
-  #This part is currently not working
-  post_observed <- posterior_emax
-  model_shapes  <- c("linear", "emax", "exponential")
-
-  # Option a) Simplified approach by using frequentist idea
-  fit_simple <- getModelFits(
-    models      = model_shapes,
-    dose_levels = dose_levels,
-posterior   = post_observed,
-    simple      = TRUE)
-
-  # Option b) Making use of the complete posterior distribution
-  fit <- getModelFits(
-    models      = model_shapes,
-    dose_levels = dose_levels,
-    posterior   = post_observed,
-    simple      = FALSE)
-
-  result_simple <- getBootstrapBands(fit_simple)
-  result <- getBootstrapBands(fit)
-  expect_type(result_simple, "list")
-  expect_type(result, "list")
-
-  result_2_simple <- getBootstrapBands(fit_simple, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3))
-  result_2 <- getBootstrapBands(fit, n_samples = 1e2, alpha = c(0.1, 0.9), avg_fit = FALSE, dose_seq = c(1, 2, 3))
-  expect_type(result_2_simple, "list")
-  expect_type(result_2, "list")
-
-  result_3_simple <- getBootstrapBands(fit_simple, dose_seq = NULL)
-  result_3 <- getBootstrapBands(fit, dose_seq = NULL)
-  expect_type(result_3_simple, "list")
-  expect_type(result_3, "list")
-})
-
diff --git a/tests/testthat/test-modelling.R b/tests/testthat/test-modelling.R
deleted file mode 100644
index e398d28..0000000
--- a/tests/testthat/test-modelling.R
+++ /dev/null
@@ -1,27 +0,0 @@
-# Test data
-test_data <- data.frame(
-  simulation = rep(1, 6),
-  dose = c(0, 1, 2, 3, 4, 5),
-  response = c(0, 1, 2, 3, 4, 5)
-)
-
-# Mock getPosterior function
-getPosterior <- function(data, prior_list, mu_hat, sd_hat) {
-  list(
-    means = c(0, 1, 2, 3, 4, 5),
-    vars = c(1, 1, 1, 1, 1, 1),
-    weights = c(1, 1, 1, 1, 1, 1)
-  )
-}
-
-# Test predictModelFit function
-test_that("predictModelFit works correctly", {
-  model_fit <- list(
-    model = "emax",
-    coeffs = c(e0 = 0, eMax = 1, ed50 = 2),
-    dose_levels = c(0, 1, 2, 3, 4, 5)
-  )
-  
-  pred_vals <- predictModelFit(model_fit)
-  expect_type(pred_vals, "double")
-})
diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R
deleted file mode 100644
index 732580d..0000000
--- a/tests/testthat/test-plot.R
+++ /dev/null
@@ -1,171 +0,0 @@
-test_that("Test plot.modelFits with different model_fits input", {
-  library(BayesianMCPMod)
-  library(clinDR)
-  library(dplyr)
-  data("metaData")
-  testdata <- as.data.frame(metaData)
-  dataset <- filter(testdata, bname == "VIAGRA")
-  histcontrol <- filter(dataset, dose == 0, primtime == 12, indication == "ERECTILE DYSFUNCTION")
-  
-  ## Create MAP Prior
-  hist_data <- data.frame(
-    trial = c(1, 2, 3, 4),
-    est = histcontrol$rslt,
-    se = histcontrol$se,
-    sd = histcontrol$sd,
-    n = histcontrol$sampsize
-  )
-  
-  sd_tot <- with(hist_data, sum(sd * n) / sum(n))
-  
-  
-  gmap <- RBesT::gMAP(
-    formula    = cbind(est, se) ~ 1 | trial,
-    weights    = hist_data$n,
-    data       = hist_data,
-    family     = gaussian,
-    beta.prior = cbind(0, 100 * sd_tot),
-    tau.dist   = "HalfNormal",
-    tau.prior  = cbind(0, sd_tot / 4)
-  )
-  
-  prior_ctr <- RBesT::robustify(
-    priormix = RBesT::automixfit(gmap),
-    weight   = 0.5,
-    mean     = 1.4,
-    sigma    = sd_tot
-  )
-  
-  # RBesT::ess(prior_ctr)
-  
-  ## derive prior for treatment
-  ## weak informative using same parameters as for robustify component
-  prior_trt <- RBesT::mixnorm(
-    comp1 = c(w = 1, m = 1.4, n = 1),
-    sigma = sd_tot,
-    param = "mn"
-  )
-  dose_levels <- c(0, 50, 100, 200)
-  ## combine priors in list
-  prior_list <- c(list(prior_ctr), rep(list(prior_trt), times = length(dose_levels[-1])))
-
-  # Pre-Specification (B)MCPMod
-  
-  ## candidate models for MCPMod
-  # linear function - no guestimates needed
-  exp <- DoseFinding::guesst(
-    d = 50,
-    p = c(0.2),
-    model = "exponential",
-    Maxd = max(dose_levels)
-  )
-  emax <- DoseFinding::guesst(
-    d = 100,
-    p = c(0.9),
-    model = "emax"
-  )
-  
-  
-  mods <- DoseFinding::Mods(
-    linear      = NULL,
-    emax        = emax,
-    exponential = exp,
-    doses       = dose_levels,
-    maxEff      = 10,
-    placEff     = 1.4
-  )
-
-  # Simulation of new trial
-  ## Note: This part will be simplified and direct results from one trial will be used
-  mods_sim <- DoseFinding::Mods(
-    emax        = emax,
-    doses       = dose_levels,
-    maxEff      = 12,
-    placEff     = 1.4
-  )
-  
-  n_patients <- c(50, 50, 50, 50)
-  data <- simulateData(
-    n_patients  = n_patients,
-    dose_levels = dose_levels,
-    sd          = sd_tot,
-    mods        = mods_sim,
-    n_sim       = 1
-  )
-  
-  data_emax <- data[, c("simulation", "dose", "emax")]
-  names(data_emax)[3] <- "response"
-  
-  posterior_emax <- getPosterior(
-    data       = data_emax,
-    prior_list = prior_list
-  )
-
-  # Evaluation of Bayesian MCPMod
-  
-  contr_mat <- DoseFinding::optContr(
-    models = mods,
-    doses  = dose_levels,
-    w      = n_patients
-  )
-  ## Calculation of critical value can be done with critVal function
-  crit_val_equal <- DoseFinding:::critVal(contr_mat$corMat, alpha = 0.05, df = 0, alternative = "one.sided")
-  crit_pval <- stats::pnorm(crit_val_equal)
-  
-  ess_prior <- round(unlist(lapply(prior_list, RBesT::ess)))
-  
-  ### Evaluation of Bayesian MCPMod
-  contr_mat_prior <- DoseFinding::optContr(
-    models = mods,
-    doses  = dose_levels,
-    w      = n_patients + ess_prior
-  )
-  
-  BMCP_result <- performBayesianMCP(
-    posteriors_list = list(posterior_emax),
-    contr_mat       = contr_mat_prior,
-    crit_prob       = crit_pval
-  )
-
-  # Model fit
-  # This part is currently not working
-  post_observed <- posterior_emax
-  model_shapes <- c("linear", "emax", "exponential")
-  
-  
-  # Option a) Simplified approach by using frequentist idea
-  fit_simple <- getModelFits(
-    models      = model_shapes,
-    dose_levels = dose_levels,
-    posterior   = post_observed,
-    simple      = TRUE
-  )
-  
-  # Option b) Making use of the complete posterior distribution
-  fit <- getModelFits(
-    models      = model_shapes,
-    dose_levels = dose_levels,
-    posterior   = post_observed,
-    simple      = FALSE
-  )
-  
-  # Test with default parameters and more models
-  plot1 <- plot.modelFits(fit)
-  expect_s3_class(plot1, "ggplot")
-  
-  # Test with cr_intv = TRUE and more models
-  plot2 <- plot.modelFits(fit, cr_intv = TRUE)
-  expect_s3_class(plot2, "ggplot")
-  
-  # Test with gAIC = FALSE and more models
-  plot3 <- plot.modelFits(fit, gAIC = FALSE)
-  expect_s3_class(plot3, "ggplot")
-  
-  # Test with avg_fit = FALSE and more models
-  plot4 <- plot.modelFits(fit, avg_fit = FALSE)
-  expect_s3_class(plot4, "ggplot")
-  
-  # Test with all non-default parameters and more models
-  plot5 <- plot.modelFits(fit, cr_intv = TRUE, gAIC = FALSE, avg_fit = FALSE)
-  expect_s3_class(plot5, "ggplot")
-})
\ No newline at end of file
diff --git a/tests/testthat/test-posterior.R b/tests/testthat/test-posterior.R
deleted file mode 100644
index 47eeb69..0000000
--- a/tests/testthat/test-posterior.R
+++ /dev/null
@@ -1,43 +0,0 @@
-test_that("getPosterior works correctly", {
-  # Prepare test data and parameters
-  data <- data.frame(simulation = rep(1, 4),
-                     dose = c(0, 1, 2, 3),
-                     response = c(10, 20, 30, 40))
-  prior_list <- list(1, 2, 3, 4)
-  mu_hat <- c(10, 20, 30, 40)
-  sd_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1)
-  
-  # Test getPosterior function
-  posterior_list <- getPosterior(data, prior_list, mu_hat, sd_hat)
-  expect_type(posterior_list, "character")
-  expect_s3_class(posterior_list, "postList")
-})
-
-test_that("getPosteriorI works correctly", {
-  # Prepare test data and parameters
-  data_i <- data.frame(dose = c(0, 1, 2, 3),
-                       response = c(10, 20, 30, 40))
-  prior_list <- list(1, 2, 3, 4)
-  mu_hat <- c(10, 20, 30, 40)
-  sd_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1)
-  
-  # Test getPosteriorI function
-  post_list <- getPosteriorI(data_i, prior_list, mu_hat, sd_hat)
-  expect_type(post_list, "character")
-  expect_s3_class(post_list, "postList")
-})
-
-test_that("summary.postList works correctly", {
-  # Prepare test data
-  post_list <- list(
-    Ctr = matrix(c(0.25, 10, 1), nrow = 3, ncol = 1),
-    DG_1 = matrix(c(0.25, 20, 2), nrow = 3, ncol = 1),
-    DG_2 = matrix(c(0.25, 30, 3), nrow = 3, ncol = 1),
-    DG_3 = matrix(c(0.25, 40, 4), nrow = 3, ncol = 1)
-  )
-  class(post_list) <- "postList"
-  
-  # Test summary.postList function
-  summary_tab <- summary.postList(post_list)
-  expect_type(summary_tab, "character")
-})
diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd
index 8e43208..8f6c247 100644
--- a/vignettes/Simulation_Example.Rmd
+++ b/vignettes/Simulation_Example.Rmd
@@ -26,15 +26,15 @@ set.seed(7015)
 # Background and data
 
 In this vignette we will show the use of the Bayesian MCPMod package for trial planning for continuous distributed data. 
-As in [link other vignette] we focus on the indication MDD and make use of historical data that is included in the clinDR package. 
+As in [the analysis example vignette](analysis_normal.html) we focus on the indication MDD and make use of historical data that is included in the clinDR package. 
 More specifically trial results for BRINTELLIX will be utilized to establish an informative prior for the control group.
 
 # Calculation of a MAP prior
-In a first step a meta analytic prior will be calculated using historical data from 5 trials (with main endpoint CfB in MADRS score after 8 weeks).
-Please note that only information from the control group will be integrated (leading to an informative multicomponent prior for the control group), while for the active groups a non-informative prior will be specified.
+In a first step a meta analytic prior will be calculated using historical data from 5 trials (with main endpoint Change from baseline in MADRS score after 8 weeks).
+Please note that only information from the control group will be integrated (leading to an informative mixture prior for the control group), while for the active groups a non-informative prior will be specified.
 
 
-```{r Calculation of a MAP prior}
+```{r Historical Data}
 data("metaData")
 testdata    <- as.data.frame(metaData)
 dataset     <- filter(testdata, bname == "BRINTELLIX")
@@ -48,19 +48,81 @@ hist_data <- data.frame(
   sd    = histcontrol$sd,
   n     = histcontrol$sampsize)
 
+sd_tot <- with(hist_data, sum(sd * n) / sum(n))
+```
+
+We will make use of the same getPriorList function as in the [analysis example vignette](analysis_normal.html) to create a MAP prior.
+
+```{r Calculation of a MAP prior, echo = FALSE}
+getPriorList <- function (
+  
+  hist_data,
+  dose_levels,
+  dose_names    = NULL,
+  robust_weight = 0.5
+  
+) {
+  
+  sd_tot <- with(hist_data, sum(sd * n) / sum(n))
+  
+  gmap <- RBesT::gMAP(
+    formula    = cbind(est, se) ~ 1 | trial,
+    weights    = hist_data$n,
+    data       = hist_data,
+    family     = gaussian,
+    beta.prior = cbind(0, 100 * sd_tot),
+    tau.dist   = "HalfNormal",
+    tau.prior  = cbind(0, sd_tot / 4))
+  
+  prior_ctr <- RBesT::automixfit(gmap)
+  
+  if (!is.null(robust_weight)) {
+    
+    prior_ctr <- suppressMessages(RBesT::robustify(
+      priormix = prior_ctr,
+      weight   = robust_weight,
+      sigma    = sd_tot))
+    
+  }
+  
+  prior_trt <- RBesT::mixnorm(
+    comp1 = c(w = 1, m = summary(prior_ctr)[1], n = 1),
+    sigma = sd_tot,
+    param = "mn")
+  
+  prior_list <- c(list(prior_ctr),
+                  rep(x     = list(prior_trt),
+                      times = length(dose_levels[-1])))
+  
+  if (is.null(dose_names)) {
+    
+    dose_names <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1])))
+    
+  }
+  
+  names(prior_list) <- dose_names
+  
+  return (prior_list)
+  
+}
+```
+
+
+```{r Evaluation of a MAP prior}
 dose_levels <- c(0, 2.5, 5, 10, 20)
 
 prior_list <- getPriorList(
   hist_data   = hist_data,
   dose_levels = dose_levels,
-  robustify_weight = 0.5)
+  robust_weight = 0.3)
 
+getESS(prior_list)
 ```
 
 # Specification of new trial design 
 
-For the hypothetical new trial we plan with 4 active dose levels $2.5, 5, 10, 20$ and we specify a broad set of potential dose-response relationships, including a linear, an exponential, an emax and a sigEMax model.  
-Furthermore we assume a maximum effect of -3 on top of control (i.e. assuming that  active treatment can reduce the MADRS score after 8 weeks by up to 15.8) and plan a trial with 80 patients for all active groups and 60 patients for control.
+For the hypothetical new trial we plan with 4 active dose levels \eqn{2.5, 5, 10, 20} and we specify a broad set of potential dose-response relationships, including a linear, an exponential, an emax and a sigEMax model.  
+Furthermore we assume a maximum effect of -3 on top of control (i.e. assuming that active treatment can reduce the MADRS score after 8 weeks by up to 15.8) and plan a trial with 80 patients for all active groups and 60 patients for control.
 ```{r}
 #Pre-Specification (B)MCPMod 
 
@@ -73,22 +135,25 @@ exp     <- DoseFinding::guesst(d     = 5,
 emax    <- DoseFinding::guesst(d     = 2.5,
                                p     = c(0.9),
                                model = "emax")
-sigemax<-  DoseFinding::guesst(d     = c(2.5,5),
-                               p     = c(0.1,0.6),
+sigemax <- DoseFinding::guesst(d     = c(2.5, 5),
+                               p     = c(0.1, 0.6),
+                               model = "sigEmax")
+sigemax2 <- DoseFinding::guesst(d     = c(2, 4),
+                               p     = c(0.3, 0.8),
                                model = "sigEmax")
-#beta <- DoseFinding::guesst(d=5, p=0.8, model="betaMod", dMax=1, scal=1.2, Maxd=20)
+
 
 mods <- DoseFinding::Mods(
   linear      = NULL,
   emax        = emax,
   exponential = exp,
-  sigEmax     = sigemax,
+  sigEmax     = rbind(sigemax, sigemax2),
   #betaMod     = beta,
   doses       = dose_levels,
   maxEff      = -3,
   placEff     = -12.8)
 
-n_patients=c(60,80,80,80,80)
+n_patients = c(60, 80, 80, 80, 80)
 ```
 
 # Calculation of success probabilities
@@ -96,21 +161,38 @@ n_patients=c(60,80,80,80,80)
 To calculate success probabilities for the different assumed dose-response models and the specified trial design we will apply the assessDesign function.  
 
 ```{r}
-success_probabilities<-assessDesign (n_patients=n_patients,
-  mods=mods,
-  prior_list=prior_list)
+success_probabilities <- assessDesign(
+  n_patients  = n_patients,
+  mods        = mods,
+  prior_list  = prior_list,
+  sd          = sd_tot)
+
 success_probabilities
 ```
 
 As alternative design we will evaluate a design with the same overall sample size but putting more patients on the highest dose group (and control).
 
 ```{r}
-success_probabilities_uneq<-assessDesign (n_patients=c(80,60,60,60,120),
-  mods=mods,
-  prior_list=prior_list)
+success_probabilities_uneq <- assessDesign(
+  n_patients  = c(80, 60, 60, 60, 120),
+  mods        = mods,
+  prior_list  = prior_list,
+  sd          = sd_tot)
 success_probabilities_uneq
 ```
 
 For this specific trial setting the adapted allocation ratio leads to increased success probabilities under all assumed dose response relationships.
 
+Instead of specifying the assumed effects via the models it is also possible to directly specify the effects for the individual dose levels via the dr_means input. This allows e.g. also the simulation of scenarios with a prior-data conflict.
+
+```{r}
+success_probabilities <- assessDesign(
+  n_patients  = c(60, 80, 80, 80, 80),
+  mods        = mods,
+  prior_list  = prior_list,
+  sd          = sd_tot,
+  dr_means = c(-12,-14,-15,-16,-17))
+success_probabilities
+```
+
 
diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd
index 8420240..870c640 100644
--- a/vignettes/analysis_normal.Rmd
+++ b/vignettes/analysis_normal.Rmd
@@ -28,17 +28,20 @@ In this vignette we will show the use of the Bayesian MCPMod package for continu
 Hereby the focus is on the utilization of an informative prior and the BayesianMCPMod evaluation of a single trial.
 We will use data that is included in the clinDR package. 
 More specifically trial results for BRINTELLIX will be used to illustrate the specification of an informative prior and the usage of such a prior for the bayesian evaluation of a (hypothetical) new trial. 
-More information around BRINTELLIX to be added...
+BRINTELLIX is a medication used to treat major depressive disorder. Various clinical trials with different dose groups (including a control group) were conducted.
 
 
 # Calculation of a MAP prior
-In a first step a meta analytic prior will be calculated using historical data from 4 trials (with main endpoint CfB in MADRS score after 8 weeks).
-Please note that only information from the control group will be integrated (leading to an informative multicomponent prior for the control group), while for the active groups a non-informative prior will be specified.
-```{r Calculation of a MAP prior}
+In a first step a meta analytic prior will be calculated using historical data from 4 trials (with main endpoint Change from baseline in MADRS score after 8 weeks).
+Please note that only information from the control group will be integrated (leading to an informative mixture prior for the control group), while for the active groups a non-informative prior will be specified.
+```{r Historical Data for Control Arm}
 data("metaData")
-testdata    <- as.data.frame(metaData)
-dataset     <- filter(testdata, bname == "BRINTELLIX")
-histcontrol <- filter(dataset, dose == 0, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER",protid!=6)
+dataset     <- filter(as.data.frame(metaData), bname == "BRINTELLIX")
+histcontrol <- filter(dataset,
+                      dose       == 0,
+                      primtime   == 8,
+                      indication == "MAJOR DEPRESSIVE DISORDER",
+                      protid     != 5)
 
 ##Create MAP Prior
 hist_data <- data.frame(
@@ -47,100 +50,181 @@ hist_data <- data.frame(
   se    = histcontrol$se,
   sd    = histcontrol$sd,
   n     = histcontrol$sampsize)
-
+```
+Here, we suggest a function to construct a list of prior distributions for the different dose groups.
+This function is adapted to the needs of this example.
+Other applications may need a different way to construct prior distributions.
+```{r Defining MAP prior function}
+getPriorList <- function (
+  
+  hist_data,
+  dose_levels,
+  dose_names    = NULL,
+  robust_weight = 0.5
+  
+) {
+  
+  sd_tot <- with(hist_data, sum(sd * n) / sum(n))
+  
+  gmap <- RBesT::gMAP(
+    formula    = cbind(est, se) ~ 1 | trial,
+    weights    = hist_data$n,
+    data       = hist_data,
+    family     = gaussian,
+    beta.prior = cbind(0, 100 * sd_tot),
+    tau.dist   = "HalfNormal",
+    tau.prior  = cbind(0, sd_tot / 4))
+  
+  prior_ctr <- RBesT::automixfit(gmap)
+  
+  if (!is.null(robust_weight)) {
+    
+    prior_ctr <- suppressMessages(RBesT::robustify(
+      priormix = prior_ctr,
+      weight   = robust_weight,
+      sigma    = sd_tot))
+    
+  }
+  
+  prior_trt <- RBesT::mixnorm(
+    comp1 = c(w = 1, m = summary(prior_ctr)[1], n = 1),
+    sigma = sd_tot,
+    param = "mn")
+  
+  prior_list <- c(list(prior_ctr),
+                  rep(x     = list(prior_trt),
+                      times = length(dose_levels[-1])))
+  
+  if (is.null(dose_names)) {
+    
+    dose_names <- c("Ctr", paste0("DG_", seq_along(dose_levels[-1])))
+    
+  }
+  
+  names(prior_list) <- dose_names
+  
+  return (prior_list)
+  
+}
+```
+With the dose levels to be investigated, the prior distribution can be constructed.
+```{r Getting the MAP prior}
 dose_levels <- c(0, 2.5, 5, 10)
 
 prior_list <- getPriorList(
   hist_data   = hist_data,
-  dose_levels = dose_levels,robustify_weight = 0.5)
+  dose_levels = dose_levels,
+  robust_weight = 0.3)
 
+getESS(prior_list)
 ```
-
 # Specifications new trial 
-We will use the trial with ct.gov number NCT00635219 as exemplary new trial. 
-To be able to apply the bayesian MCPMod approach, candidate models need to be specified. Since there are only 3 active dose levels we will limit the set of candidate models to a linear, an exponential and an emax model.
+We will use the trial with ct.gov number NCT00735709 as exemplary new trial. 
+To be able to apply the bayesian MCPMod approach, candidate models need to be specified.
+Since there are only 3 active dose levels we will limit the set of candidate models to a linear, an exponential, and an emax model.
+Note that the linear candidate model does not require a guesstimate.
 
 ```{r Pre-Specification of candidate models}
-#Pre-Specification (B)MCPMod 
-
-## candidate models for MCPMod
-# linear function - no guestimates needed
-exp     <- DoseFinding::guesst(d     = 5,
-                               p     = c(0.2),
-                               model = "exponential",
-                               Maxd  = max(dose_levels))
-emax    <- DoseFinding::guesst(d     = 2.5,
-                               p     = c(0.9),
-                               model = "emax")
+exp_guesst  <- DoseFinding::guesst(d     = 5,
+                                   p     = c(0.2),
+                                   model = "exponential",
+                                   Maxd  = max(dose_levels))
 
+emax_guesst <- DoseFinding::guesst(d     = 2.5,
+                                   p     = c(0.9),
+                                   model = "emax")
 
 mods <- DoseFinding::Mods(
   linear      = NULL,
-  emax        = emax,
-  exponential = exp,
+  emax        = emax_guesst,
+  exponential = exp_guesst,
   doses       = dose_levels,
   maxEff      = -1,
   placEff     = -12.8)
 
-new_trial <- filter(dataset, primtime == 8, indication == "MAJOR DEPRESSIVE DISORDER",protid==6)
-n_patients <- c(150, 142, 147, 149)
+new_trial  <- filter(dataset,
+                     primtime   == 8,
+                     indication == "MAJOR DEPRESSIVE DISORDER",
+                     protid     == 5)
+
+n_patients <- c(128, 124, 129, 122)
 
 ```
 
 # Combination of prior information and trial results
 
-As outlined in citePaper, in a first step the posterior is calculated combining the prior information with the estimated results of the new trial.
+As outlined in Fleischer et al. [Bayesian MCPMod. Pharmaceutical Statistics. 2022; 21(3): 654-670.], in a first step the posterior is calculated combining the prior information with the estimated results of the new trial. Via the summary function it is possible to print out the summary information of the posterior distributions.
+
 
 ```{r Trial results}
-posterior <- getPosterior(prior=prior_list,
-                               mu_hat = new_trial$rslt,
-                               sd_hat = new_trial$se)
+
+posterior <- getPosterior(prior  = prior_list,
+                          mu_hat = new_trial$rslt,
+                          se_hat = new_trial$se,
+                          calc_ess = TRUE)
+
+summary(posterior)
+
 
 ```
 
 # Execution of Bayesian MCPMod Test step
-For the execution of the testing step of bayesian MCPMod a critical value (on the probability scale) will be determined for a given alpha level. In addition a contrast matrix is created. Please note that there are different possibilities how to generate contrasts.
-This information is then used as input for the bayesian MCP testing function. 
+For the execution of the testing step of Bayesian MCPMod a critical value (on the probability scale) will be determined for a given alpha level. This critical value is calculated using (the re-estimated) contrasts for the frequentist MCPMod to ensure that when using non-informative priors the actual error level for falsely declaring a significant trial in the Bayesian MCPMod is controlled (by the specified alpha level). 
 
-```{r Execution of Bayesian MCPMod Test step}
-crit_pval <- getCritProb(
-  mods           = mods,
-  dose_levels    = dose_levels,
-  dose_weights   = n_patients,
-  alpha_crit_val = 0.1)
-
-contr_mat_prior <- getContrMat(
-  mods           = mods,
-  dose_levels    = dose_levels,
-  dose_weights   = n_patients,
-  prior_list     = prior_list)
-
-#This would be the most reasonable output, but since it is not exported it is currently not usable.
-#BMCP_result<-BayesMCPi(posterior_i = posterior,
-     #     contr_mat       = contr_mat_prior, 
-      #    crit_prob       = crit_pval)
+A pseudo-optimal contrast matrix is generated based on the variability of the posterior distribution (see Fleischer et al. 2022 for more details). Please note that there are different possibilities how to establish the contrasts.
 
-BMCP_result <- performBayesianMCP(
-  posteriors_list = posterior,
-  contr_mat       = contr_mat_prior, 
- crit_prob       = crit_pval)
 
-#BMCP_result2 <- performBayesianMCPMod(
-#    posteriors_list = posterior_emax,
-#    contr_mat       = contr_mat_prior, 
-#    crit_prob       = crit_pval)
 
-#BMCP_result2
+```{r Preparation of input for Bayesian MCPMod Test step}
+crit_pval <- getCritProb(
+    mods           = mods,
+    dose_levels    = dose_levels,
+    se_new_trial   = new_trial$se,
+    alpha_crit_val = 0.05)
+
+contr_mat<- getContr(
+ mods           = mods,
+ dose_levels    = dose_levels,
+ sd_posterior   = summary(posterior)[,2])
+
+#Alternative ways to come up with a contrast would be to use 
+#i) the frequentist contrast 
+#contr_mat_prior <- getContr(
+# mods           = mods,
+# dose_levels    = dose_levels,
+# dose_weights   = n_patients,
+# prior_list     = prior_list)
+# ii) re-estimated frequentist contrasts
+# contr_mat_prior <- getContr(
+# mods           = mods,
+# dose_levels    = dose_levels,
+# se_new_trial   = new_trial$se)
+# iii)  Bayesian approach using number of patients for new trial and prior distribution  
+#contr_mat_prior <- getContr(
+#  mods           = mods,
+#  dose_levels    = dose_levels,
+#  dose_weights   = n_patients,
+#  prior_list     = prior_list)
+```
+The bayesian MCP testing step is then executed based on the posterior information, the provided contrasts and the (multiplicity adjusted) critical value.   
+```{r  Executionr of Bayesian MCPMod Test step}
+BMCP_result <- performBayesianMCP(
+  posterior_list = posterior,
+  contr          = contr_mat, 
+  crit_prob_adj  = crit_pval)
 
-BMCP_result[1,1:5]
+BMCP_result
 ```
-The testing step is significant indicating a non-flat dose-response shape. In detail the p-values for the emax model is the most significant one.
+
+The testing step is significant indicating a non-flat dose-response shape. In detail, all 3 models are significant and the p-value for the emax model is the most significant one.
 
 # Model fitting and visualization of results
-In the model fitting step the posterior distribution is used as basis. Both simplified and full fitting are performed.
-Please note that all models are fitted to the data even though they were not significant in the testing step.
-For the plotting bootstrap based credible intervals (for 50% and 95%) are shown as well.
 
+In the model fitting step the posterior distribution is used as basis. Both simplified and full fitting are performed. 
+For the simplified fit the multivariate normal distribution of the control group is approximated and reduced by a one-dimensional normal distribution.The actual fit (on this approximated posterior distribution) is then performed using generalized least squares criterion.
+For the full fitting the (multi-dimensional) non-linear optimization problem is adressed via the Nelder Mead algorithm. 
+
+The output of the fit includes information about the predicted effects for the included dose levels the generalized AIC (and the corresponding weights). For the considered case the simplified and the full fit are very similar.
 ```{r Model fitting}
 #Model fit
 post_observed <- posterior
@@ -159,12 +243,40 @@ fit <- getModelFits(
   dose_levels = dose_levels,
   posterior   = post_observed,
   simple      = FALSE)
+```
+
+Via the predict function one can also receive estimates for dose levels that were not included in the trial.
+
+It is possible to plot the fitted dose response models and an AIC based average model (black lines). To assess the uncertainty one can in addition visualize credible bands (orange shaded areas, the default is set to 50% and 95%). These credible bands are calculated as follows.
+Samples from the posterior distribution are drawn and for every sample the simplified fitting step and a prediction is performed. These fits are then used to identify and visualize the specified quantiles. 
+The bootstrap based quantiles can also be directly calculated and displayed via the getbootstrapQuantiles function.
+ 
+
+```{r Visualization and interpretation of the results}
+predict(fit, doses = c(0, 2.5,4, 5,7, 10))
+
 
-plot(fit,cr_bands=TRUE)
+plot(fit, cr_bands = TRUE)
 plot(fit_simple, cr_bands = TRUE)
 
+
+getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 2.5,4, 5,7, 10))
 ```
+
+Technical note: The median quantile of the bootstrap based procedure is not necessary similar to the main model fit, as they are derived via different procedures (main fit, i.e. black lines in the plot, show the best fit of a certain model based on minimizing the residuals for the posterior distribution while bootstrap based 50% quantile indicates the median fit of the random sampling and fitting procedure). 
+
 # Additional notes
-TBD, whether certain wrapper functions should be mentioned.
+It is also possible to perform the testing and modelling step in a combined fashion via the performBayesianMCPMod function.
+
+```{r}
+overall_result<-performBayesianMCPMod(
+      posterior_list   = posterior,
+      contr            = contr_mat, 
+      crit_prob_adj    = crit_pval,
+      simple           = FALSE)
+
+overall_result$Mod
+```
+