diff --git a/R/modelling.R b/R/modelling.R
index 22a33e6..2e66e4e 100644
--- a/R/modelling.R
+++ b/R/modelling.R
@@ -49,7 +49,7 @@ getModelFitSimple <- function (
     resp  = summary.postList(posterior)[, 1],
     S     = diag(summary.postList(posterior)[, 2]^2),
     model = model,
-    type  = doFit_ID$VALUES$FIT_TYPE,
+    type  = "general",
     bnds  = DoseFinding::defBnds(mD = max(dose_levels))[[model]])
   
   pred_vals  <- stats::predict(fit, predType = "ls-means")
@@ -190,31 +190,37 @@ predictModelFit <- function (
   
   pred_vals <- switch (
     model_fit$model,
-    "emax" = {DoseFinding::emax(doses,
-                                model_fit$coeffs["e0"],
-                                model_fit$coeffs["eMax"],
-                                model_fit$coeffs["ed50"])},
-    "sigEmax" = {DoseFinding::sigEmax(doses,
-                                      model_fit$coeffs["e0"],
-                                      model_fit$coeffs["eMax"],
-                                      model_fit$coeffs["ed50"],
-                                      model_fit$coeffs["h"])},
-    "exponential" = {DoseFinding::exponential(doses,
-                                              model_fit$coeffs["e0"],
-                                              model_fit$coeffs["e1"],
-                                              model_fit$coeffs["delta"])},
-    "quadratic" = {DoseFinding::quadratic(doses,
-                                          model_fit$coeffs["e0"],
-                                          model_fit$coeffs["b1"],
-                                          model_fit$coeffs["b2"])},
-    "linear" = {DoseFinding::linear(doses,
-                                    model_fit$coeffs["e0"],
-                                    model_fit$coeffs["delta"])},
-    "logistic" = {DoseFinding::logistic(doses,
-                                        model_fit$coeffs["e0"],
-                                        model_fit$coeffs["eMax"],
-                                        model_fit$coeffs["ed50"],
-                                        model_fit$coeffs["delta"])},
+    "emax" = {DoseFinding::emax(
+      doses,
+      model_fit$coeffs["e0"],
+      model_fit$coeffs["eMax"],
+      model_fit$coeffs["ed50"])},
+    "sigEmax" = {DoseFinding::sigEmax(
+      doses,
+      model_fit$coeffs["e0"],
+      model_fit$coeffs["eMax"],
+      model_fit$coeffs["ed50"],
+      model_fit$coeffs["h"])},
+    "exponential" = {DoseFinding::exponential(
+      doses,
+      model_fit$coeffs["e0"],
+      model_fit$coeffs["e1"],
+      model_fit$coeffs["delta"])},
+    "quadratic" = {DoseFinding::quadratic(
+      doses,
+      model_fit$coeffs["e0"],
+      model_fit$coeffs["b1"],
+      model_fit$coeffs["b2"])},
+    "linear" = {DoseFinding::linear(
+      doses,
+      model_fit$coeffs["e0"],
+      model_fit$coeffs["delta"])},
+    "logistic" = {DoseFinding::logistic(
+      doses,
+      model_fit$coeffs["e0"],
+      model_fit$coeffs["eMax"],
+      model_fit$coeffs["ed50"],
+      model_fit$coeffs["delta"])},
     {stop(GENERAL$ERROR$MODEL_OPTIONS)})
   
   return (pred_vals)
diff --git a/R/plot.R b/R/plot.R
index d15ca48..64d675c 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -1,53 +1,111 @@
-plot.estMod <- function (
+plot.modelFits <- function (
     
-  est_mod
-  # dose_levels
-  # posteriors = posterior_linear[[2]]
+  model_fits,
+  CrI     = FALSE,
+  gAIC    = TRUE,
+  avg_fit = TRUE
   
 ) {
   
-  model <- est_mod$model
-  theta <- est_mod$fit$solution
-  
-  dose  <- seq(min(dose_levels), max(dose_levels), length.out = 1e3)
-  
-  switch(model,
-         "emax" = {
-           resp_expr <- quote(theta[1] + (theta[2] * dose) / (theta[3] + dose))},
-         "sigEmax" = {
-           resp_expr <- quote(theta[1] + (theta[2] * dose^theta[4]) / (theta[3]^theta[4] + dose^theta[4]))},
-         "exponential" = {
-           resp_expr <- quote(theta[1] + theta[2] * (exp(dose / theta[3]) - 1))},
-         "quadratic" = {
-           resp_expr <- quote(theta[1] + theta[2] * dose + theta[3] * dose^2)},
-         "linear" = {
-           resp_expr  <- quote(theta[1] + theta[2] * dose)},
-         "logistic" = {
-           resp_expr <- quote(theta[1] + theta[2] / (1 + exp((theta[3] - dose) / theta[4])))},
-         {
-           stop(GENERAL$ERROR$MODEL_OPTIONS)}
-  )
-  
-  df <- data.frame(dose = dose, response = eval(resp_expr))
-  
-  data.frame(dose = dose_levels, obs = )
-  
-  plt <- ggplot2::ggplot(data = df) +
-    ggplot2::geom_line(ggplot2::aes(dose, response)) +
-    ggplot2::geom_point()
-  
-  return(plt)
-  
-}
-
-plot.estMods <- function (
+  plot_resolution <- 1e3
+  
+  dose_levels  <- model_fits[[1]]$dose_levels
+  post_summary <- summary.postList(attr(model_fits, "posterior"))
+  doses        <- seq(from = min(dose_levels),
+                      to   = max(dose_levels), length.out = plot_resolution)
+  
+  preds_models <- sapply(model_fits, predictModelFit, doses = doses)
+  model_names  <- names(model_fits)
+  
+  if (avg_fit) {
+    
+    mod_weigts <- sapply(model_fits, function (x) x$model_weight)
+    avg_mod    <- preds_models %*% mod_weigts
+    
+    preds_models <- cbind(preds_models, avg_mod)
+    model_names  <- c(model_names, "averageModel")
     
-  est_mods
+  }
   
-) {
+  gg_data      <- data.frame(
+    dose_levels = rep(doses, length(model_names)),
+    fits        = as.vector(preds_models),
+    models      = rep(factor(model_names,
+                             levels = c("linear", "emax", "exponential",
+                                        "sigEmax", "logistic", "quadratic",
+                                        "averageModel")),
+                      each = plot_resolution))
   
-  plts <- lapply(est_mods, plot)
+  if (gAIC) {
+    
+    g_AICs     <- sapply(model_fits, function (x) x$gAIC)
+    label_gAUC <- paste("AIC:", round(g_AICs, digits = 1))
+    
+    if (avg_fit) {
+      
+      mod_weigts  <- sort(mod_weigts, decreasing = TRUE)
+      paste_names <- names(mod_weigts) |>
+        gsub("exponential", "exp", x = _) |>
+        gsub("quadratic",  "quad", x = _) |>
+        gsub("linear",      "lin", x = _) |>
+        gsub("logistic",    "log", x = _) |>
+        gsub("sigEmax",    "sigE", x = _)
+      
+      label_avg  <- paste0(paste_names, "=", round(mod_weigts, 1),
+                           collapse = ", ")
+      label_gAUC <- c(label_gAUC, label_avg)
+      
+    }
+    
+  }
   
+  plts <- ggplot2::ggplot() + 
+    ## Layout etc.
+    ggplot2::theme_bw() + 
+    ggplot2::labs(x = "Dose",
+                  y = "Model Fits") +
+    ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
+                   panel.grid.minor = ggplot2::element_blank()) +
+    ## gAIC
+    {if (gAIC) {
+      
+      ggplot2::geom_text(
+        data    = data.frame(
+          models = unique(gg_data$models),
+          label  = label_gAUC),
+        mapping = ggplot2::aes(label = label_gAUC),
+        x = -Inf, y = Inf, hjust = "inward", vjust = "inward",
+        size = 3)
+      
+      }
+    } + 
+    ## Posterior Credible Intervals
+    {if (CrI) {
+      
+      ggplot2::geom_errorbar(
+        data    = data.frame(x    = dose_levels,
+                             ymin = post_summary[, 3],
+                             ymax = post_summary[, 5]),
+        mapping = ggplot2::aes(x    = x,
+                               ymin = ymin,
+                               ymax = ymax),
+        width = 0, alpha = 0.5)
+      
+      }
+    } + 
+    ## Posterior Medians
+    ggplot2::geom_point(
+      data    = data.frame(dose_levels = dose_levels,
+                           fits        = post_summary[, 4]),
+      mapping = ggplot2::aes(dose_levels, fits),
+      size    = 2) +
+    ## Fitted Models
+    ggplot2::geom_line(
+      data    = gg_data,
+      mapping = ggplot2::aes(dose_levels, fits)) + 
+    ## Faceting
+    ggplot2::facet_wrap(~ models)
   
+  return (plts)
   
-}
+}
\ No newline at end of file
diff --git a/R/posterior.R b/R/posterior.R
index a47e1ff..42d1026 100644
--- a/R/posterior.R
+++ b/R/posterior.R
@@ -7,26 +7,56 @@
 getPosterior <- function(
     
   data,
-  prior_list
+  prior_list,
+  mu_hat = NULL,
+  sd_hat = NULL
   
 ) {
   
-  lapply(split(data, data$simulation), getPosteriorI, prior_list = prior_list)
+  posterior_list <- lapply(split(data, data$simulation), getPosteriorI,
+                           prior_list = prior_list,
+                           mu_hat     = mu_hat,
+                           sd_hat     = sd_hat)
+  
+  if (length(posterior_list) == 1) {
+    
+    posterior_list <- posterior_list[[1]]
+    
+  }
+  
+  return (posterior_list)
   
 }
 
 getPosteriorI <- function(
     
   data_i,
-  prior_list
+  prior_list,
+  mu_hat = NULL,
+  sd_hat = NULL
   
 ) {
   
-  anova_res  <- lm(data_i$response ~ factor(data_i$dose) - 1) # take mean & var out in separate function
-  anova_mean <- summary(anova_res)$coefficients[, 1]          #
-  anova_se   <- summary(anova_res)$coefficients[, 2]          #
+  if (is.null(mu_hat) && is.null(sd_hat)) {
+    
+    anova_res <- lm(data_i$response ~ factor(data_i$dose) - 1)
+    mu_hat    <- summary(anova_res)$coefficients[, 1]
+    sd_hat    <- summary(anova_res)$coefficients[, 2]
+    
+  } else if (!is.null(mu_hat) && !is.null(sd_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) == nrow(sd_hat))
+    
+  } else {
+    
+    stop ("Both mu_hat and S_hat must be provided.")
+    
+  }
   
-  post_list <- mapply(RBesT::postmix, prior_list, m = anova_mean, se = anova_se)
+  post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = sd_hat)
   
   names(post_list) <- c("Ctr", paste0("DG_", seq_along(post_list[-1])))
   class(post_list) <- "postList"
@@ -55,11 +85,12 @@ getPostCombsI <- function (
 
 summary.postList <- function (
     
-  post_list
+  post_list,
+  ...
   
 ) {
   
-  summary_list        <- lapply(post_list, summary)
+  summary_list        <- lapply(post_list, summary, ...)
   names(summary_list) <- names(post_list)
   summary_tab         <- do.call(rbind, summary_list)
   
diff --git a/vignettes/.gitignore b/vignettes/.gitignore
new file mode 100644
index 0000000..097b241
--- /dev/null
+++ b/vignettes/.gitignore
@@ -0,0 +1,2 @@
+*.html
+*.R
diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd
new file mode 100644
index 0000000..57a74f6
--- /dev/null
+++ b/vignettes/analysis_normal.Rmd
@@ -0,0 +1,176 @@
+---
+title: "Analysis Example of Bayesian MCPMod for Continuous Data"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Analysis Example of Bayesian MCPMod for Continuous Data}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(BayesianMCPMod)
+library(clinDR)
+library(dplyr)
+```
+This vignette shows how to apply the BayesianMCPMod package.
+
+# i) Calculation of a MAP prior
+In a first step a meta analytic prior will be calculated using historical data.
+```{r}
+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])))
+```
+# ii) Pre-Specification of candidate models 
+```{r}
+#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)
+```
+# iii) Trial results
+```{r}
+#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)
+
+```
+
+# iv) Execution of Bayesian MCPMod Test step
+```{r}
+#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 <- BMCPMod(
+  posteriors_list = list(posterior_emax),
+  contr_mat       = contr_mat_prior, 
+  crit_prob       = crit_pval)
+
+BMCP_result
+```
+
+# v) Model fitting and Plotting
+In the model fitting step the posterior distribution is used as basis. 
+```{r}
+#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)
+```
+
+```{r}
+plot.modelFits(fit,CrI= TRUE)
+plot.modelFits(fit_simple, CrI= TRUE)
+```
+For the plotting the credible intervals are shown as well.