diff --git a/DESCRIPTION b/DESCRIPTION
index e15756a..ec687a7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: BayesianMCPMod
 Title: Bayesian MCPMod
-Version: 0.1.3-3
+Version: 0.1.4-3
 Authors@R: c(
     person("Boehringer Ingelheim Pharma GmbH & Co. KG",
            role = c("cph", "fnd")),
@@ -11,8 +11,8 @@ Authors@R: c(
            role = "aut",
            email = "lars.andersen@boehringer-ingelheim.com"),
     person("Steven", "Brooks",
-           role = "aut",
-           email = "lars.andersen@boehringer-ingelheim.com"),
+           role = "ctb",
+           email = "steven.brooks@boehringer-ingelheim.com"),
     person("Sebastian", "Bossert",
            role = "aut",
            email = "sebastian.bossert@boehringer-ingelheim.com")
@@ -32,7 +32,8 @@ Imports:
     clinDR,
     knitr,
     rmarkdown,
-    dplyr
+    dplyr,
+    checkmate
 Suggests:
     testthat (>= 3.0.0)
 Depends:
diff --git a/R/BMCPMod.R b/R/BMCPMod.R
index 4dafc4e..39f6e00 100644
--- a/R/BMCPMod.R
+++ b/R/BMCPMod.R
@@ -24,7 +24,6 @@
 #'                    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,
@@ -53,6 +52,15 @@ assessDesign <- function (
   
 ) {
   
+  checkmate::assert_vector(n_patients, len = length(attr(mods, "doses")), any.missing = FALSE)
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
+  # sensitive to how DoseFinding labels their attributes for "Mods" class
+  checkmate::check_double(n_sim, lower = 1, upper = Inf)
+  checkmate::check_double(alpha_crit_val, lower = 0, upper = 1)
+  checkmate::check_logical(simple)
+  # TODO: check that prior_list has 'sd_tot' attribute, and that it's numeric
+  
   dose_levels <- attr(mods, "doses")
   
   data <- simulateData(
@@ -172,6 +180,11 @@ getContr <- function (
   
 ) {
   
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(attr(mods, "doses")))
+  checkmate::check_double(dose_weights, any.missing = FALSE, len = length(attr(mods, "doses")))
+  checkmate::check_list(prior_list, names = "named", len = length(attr(mods, "doses")), any.missing = FALSE)
+  
   # frequentist & re-estimation
   if (!is.null(se_new_trial) & 
       is.null(dose_weights) & is.null(prior_list) & is.null(sd_posterior)) {
@@ -267,6 +280,11 @@ getCritProb <- function (
   
 ) {
   
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(dose_weights))
+  checkmate::check_double(dose_weights, any.missing = FALSE, len = length(dose_levels))
+  checkmate::check_double(alpha_crit_val, lower = 0, upper = 1)
+  
   contr <- getContr(mods           = mods,
                     dose_levels    = dose_levels ,
                     dose_weights   = dose_weights,
@@ -329,6 +347,17 @@ performBayesianMCPMod <- function (
   
 ) {
   
+  checkmate::check_class(posterior_list, "postList")
+  checkmate::check_class(contr, "optContr")
+  checkmate::check_class(crit_prob_adj, "numeric")
+  checkmate::check_logical(simple)
+  
+  if (inherits(posterior_list,  "postList")) {
+    
+    posterior_list <- list(posterior_list)
+    
+  }
+  
   if (inherits(posterior_list,  "postList")) {
     
     posterior_list <- list(posterior_list)
@@ -453,6 +482,12 @@ performBayesianMCP <- function(
   crit_prob_adj
   
 ) {
+  
+  checkmate::check_class(posterior_list, "postList")
+  checkmate::check_class(contr, "optContr")
+  checkmate::check_class(crit_prob_adj, "numeric")
+  checkmate::check_numeric(crit_prob_adj, lower = 0, upper = Inf)
+  
   if (inherits(posterior_list,  "postList")) {
     
     posterior_list <- list(posterior_list)
diff --git a/R/modelling.R b/R/modelling.R
index ed50b13..43eab38 100644
--- a/R/modelling.R
+++ b/R/modelling.R
@@ -45,6 +45,10 @@ getModelFits <- function (
   
 ) {
   
+  checkmate::check_list(models, any.missing = FALSE)
+  checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(models))
+  checkmate::check_class(posterior, "postList")
+  checkmate::check_logical(simple)
   models      <- unique(gsub("\\d", "", models))
   
   getModelFit <- ifelse(simple, getModelFitSimple, getModelFitOpt)
@@ -137,7 +141,7 @@ getModelFitOpt <- function (
               ub     <- c(Inf, Inf, 1.5 * max(dose_levels), 0.5 * max(dose_levels))
               expr_i <- quote(sum((post_combs$means[i, ] - (theta[1] + theta[2] / (1 + exp((theta[3] - dose_levels) / theta[4]))))^2 / (post_combs$vars[i, ])))},
             {
-              stop (GENERAL$ERROR$MODEL_OPTIONS)})
+              stop ("error")})
     
     simple_fit <- getModelFitSimple(
       model       = model,
@@ -247,7 +251,7 @@ predictModelFit <- function (
       model_fit$coeffs["eMax"],
       model_fit$coeffs["ed50"],
       model_fit$coeffs["delta"])},
-    {stop(GENERAL$ERROR$MODEL_OPTIONS)})
+    {stop("error")})
   
   return (pred_vals)
   
diff --git a/R/plot.R b/R/plot.R
index b53821c..a077b71 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -43,6 +43,15 @@ plot.modelFits <- function (
   ...
   
 ) {
+  checkmate::check_class(x, "modelFits")
+  checkmate::check_logical(gAIC)
+  checkmate::check_logical(avg_fit)
+  checkmate::check_logical(cr_intv)
+  checkmate::check_double(alpha_CrI, lower = 0, upper = 1)
+  checkmate::check_logical(cr_bands)
+  checkmate::check_double(alpha_CrB, lower = 0, upper = 1, len = 2)
+  checkmate::check_integer(n_bs_smpl, lower = 1, upper = Inf)
+  checkmate::check_string(acc_color, na.ok = TRUE)
   
   plot_res   <- 1e2
   model_fits <- x
diff --git a/R/posterior.R b/R/posterior.R
index 6693fbb..1aaa41d 100644
--- a/R/posterior.R
+++ b/R/posterior.R
@@ -1,4 +1,3 @@
-
 #' @title getPosterior
 #' 
 #' @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.
@@ -36,6 +35,12 @@ getPosterior <- function(
   calc_ess = FALSE
   
 ) {
+  checkmate::check_data_frame(data, null.ok = TRUE)
+  checkmate::check_list(prior_list, names = "named", any.missing = FALSE)
+  checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE)
+  checkmate::check_double(mu_hat, null.ok = TRUE, lower = -Inf, upper = Inf)
+  checkmate::check_vector(se_hat, any.missing = FALSE, null.ok = TRUE)
+  checkmate::check_double(se_hat, null.ok = TRUE, lower = 0, upper = Inf)
   
   if (!is.null(mu_hat) && !is.null(se_hat) && is.null(data)) {
     
@@ -55,7 +60,7 @@ getPosterior <- function(
     stop ("Either 'data' or 'mu_hat' and 'se_hat' must not be NULL.")
     
   }
- 
+  
   if (length(posterior_list) == 1) {
     
     posterior_list <- posterior_list[[1]]
@@ -75,9 +80,19 @@ getPosteriorI <- function(
   calc_ess = FALSE
   
 ) {
+
+  checkmate::check_data_frame(data_i, null.ok = TRUE)
+  checkmate::check_list(prior_list, names = "named", any.missing = FALSE)
+  checkmate::check_vector(mu_hat, any.missing = FALSE, null.ok = TRUE)
+  checkmate::check_double(mu_hat, null.ok = TRUE, lower = -Inf, upper = Inf)
+  checkmate::check_vector(se_hat, any.missing = FALSE, null.ok = TRUE)
+  checkmate::check_double(se_hat, null.ok = TRUE, lower = 0, upper = Inf)
   
   if (is.null(mu_hat) && is.null(se_hat)) {
-    
+    checkmate::check_data_frame(data_i, null.ok = FALSE)
+    checkmate::assert_names(names(data_i), must.include = "response")
+    checkmate::assert_names(names(data_i), must.include = "dose")
+
     anova_res <- stats::lm(data_i$response ~ factor(data_i$dose) - 1)
     mu_hat    <- summary(anova_res)$coefficients[, 1]
     se_hat    <- summary(anova_res)$coefficients[, 2]
@@ -92,7 +107,6 @@ getPosteriorI <- function(
   } else {
     
     stop ("Both mu_hat and se_hat must be provided.")
-    
   }
   
   post_list <- mapply(RBesT::postmix, prior_list, m = mu_hat, se = se_hat,
diff --git a/R/simulation.R b/R/simulation.R
index c4bd480..50c24dc 100644
--- a/R/simulation.R
+++ b/R/simulation.R
@@ -48,6 +48,13 @@ simulateData <- function(
   
 ) {
   
+  checkmate::check_vector(n_patients, any.missing = FALSE, len = length(dose_levels))
+  checkmate::check_double(dose_levels, lower = 0, any.missing = FALSE, len = length(n_patients))
+  checkmate::check_double(sd, len = 1, null.ok = FALSE, lower = 0, upper = Inf)
+  checkmate::check_class(mods, classes = "Mods")
+  checkmate::check_numeric(n_sim, lower = 0, upper = Inf, len = 1)
+  checkmate::check_string(true_model, null.ok = TRUE)
+  
   if (!is.null(true_model)) {
     
     n_sim            <- 1
@@ -106,7 +113,7 @@ getModelData <- function (
   
 ) {
   
-  model_data              <- sim_data[, c("simulation", "dose", model_name)]
+  model_data <- sim_data[, c("simulation", "dose", model_name)]
   colnames(model_data)[3] <- "response"
   
   return (model_data)
diff --git a/man/assessDesign.Rd b/man/assessDesign.Rd
index 1f0a539..9696f54 100644
--- a/man/assessDesign.Rd
+++ b/man/assessDesign.Rd
@@ -54,7 +54,6 @@ prior_list<-list(Ctrl   = RBesT::mixnorm(comp1 = c(w = 1, m = 0, s = 12), sigma
                    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,
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
new file mode 100644
index 0000000..2a7bda7
--- /dev/null
+++ b/tests/testthat/setup.R
@@ -0,0 +1,156 @@
+#' @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 robust_weight needs to be provided as a numeric
+#' value for the weight of the robustification component
+#'
+getPriorList <- function (
+    
+  hist_data,
+  dose_levels,
+  dose_names       = NULL,
+  robust_weight
+  
+) {
+  
+  checkmate::check_data_frame(hist_data)
+  checkmate::assert_double(dose_levels, lower = 0, any.missing = FALSE)
+  checkmate::check_string(dose_names, null.ok = TRUE)
+  checkmate::check_vector(dose_names, null.ok = TRUE, len = length(dose_levels))
+  checkmate::check_numeric(robust_weight, len = 1, null.ok = FALSE)
+  
+  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)
+  
+  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)
+  
+}
+
+
+getPostProb <- function (
+    
+  contr_j,     # j: dose level
+  post_combs_i # i: simulation outcome
+  
+) {
+  
+  ## Test statistic = sum over all components of
+  ## posterior weight * normal probability distribution of
+  ## critical values for doses * estimated mean / sqrt(product of critical values for doses)
+  
+  ## Calculation for each component of the posterior
+  contr_theta   <- apply(post_combs_i$means, 1, `%*%`, contr_j)
+  contr_var     <- apply(post_combs_i$vars, 1, `%*%`, contr_j^2)
+  contr_weights <- post_combs_i$weights
+  
+  ## P(c_m * theta > 0 | Y = y) for a shape m (and dose j)
+  post_probs <- sum(contr_weights * stats::pnorm(contr_theta / sqrt(contr_var)))
+  
+  return (post_probs)
+  
+}
+
+# Create minimal test case
+n_hist_trials = 2
+
+hist_data <- data.frame(
+  trial = seq(1, n_hist_trials, 1),
+  est   = rep(1, n_hist_trials),
+  se    = rep(1, n_hist_trials),
+  sd    = rep(1, n_hist_trials),
+  n     = rep(1, n_hist_trials)
+)
+
+n_patients <- c(2, 1)
+dose_levels <- c(0, 2.5)
+mean <- c(8, 12)
+sd <- c(0.5, 0.8)
+
+mods <- DoseFinding::Mods(
+  linear = NULL, 
+  doses = dose_levels
+)
+
+
+prior_list <- getPriorList(
+  hist_data   = hist_data,
+  dose_levels = dose_levels,
+  robust_weight = 0.5
+)
+
+n_sim = 1
+alpha_crit_val = 0.05
+simple = TRUE
+
+data <- simulateData(
+  n_patients  = n_patients,
+  dose_levels = dose_levels,
+  sd          = sd,
+  mods        = mods,
+  n_sim       = n_sim
+)
+
+posterior_list <- getPosterior(
+  data = getModelData(data, names(mods)[1]),
+  prior_list = prior_list
+)
+
+contr_mat = getContr(
+  mods = mods, 
+  dose_levels = dose_levels, 
+  dose_weights = n_patients,
+  prior_list = prior_list
+)
+
+crit_pval = getCritProb(
+  mods = mods, 
+  dose_levels = dose_levels, 
+  dose_weights = n_patients, 
+  alpha_crit_val = alpha_crit_val
+)
+
+# eval_design <- assessDesign(
+#   n_patients = n_patients, 
+#   mods = mods, 
+#   prior_list = prior_list,
+#   n_sim = n_sim,
+#   alpha_crit_val = alpha_crit_val,
+#   simple = TRUE
+# )
diff --git a/tests/testthat/test-BMCPMod.R b/tests/testthat/test-BMCPMod.R
new file mode 100644
index 0000000..5ad82bf
--- /dev/null
+++ b/tests/testthat/test-BMCPMod.R
@@ -0,0 +1,296 @@
+##########################
+# Tests for assessDesign #
+##########################
+
+test_that("base case input throws no error and has correct properties", {
+  
+  expect_no_error(
+    eval_design <- assessDesign(
+      n_patients = n_patients, 
+      mods = mods, 
+      sd = sd,
+      prior_list = prior_list,
+      n_sim = n_sim,
+      alpha_crit_val = alpha_crit_val,
+      simple = TRUE
+    )
+  )
+  
+  # assessDesign should give results for each model in mods
+  expect_equal(
+    names(eval_design), names(mods)
+  )
+  
+  # assessDesign result should have rows = n_sim
+  expect_equal(
+    attr(eval_design$linear$BayesianMCP, "dim")[1],
+    n_sim
+  )
+  
+  # assessDesign result (in this base case) should have crit_prob = 1 - alpha_crit_val
+  expect_equal(
+    attr(eval_design$linear$BayesianMCP, "crit_prob"),
+    1 - alpha_crit_val
+  )
+  
+})
+
+
+### n_patients param ###
+
+test_that("assessDesign validates n_patients parameter input and give appropriate error messages", {
+  
+  # assertions that aren't tested here for sake of brevity
+    # n_patients should be a non-NULL numeric vector
+  
+  expect_error(
+    assessDesign(n_patients = n_patients[-1], sd = sd, mods = mods, prior_list = prior_list, n_sim = n_sim)
+  )
+  
+  expect_error(
+    assessDesign(n_patients = rep(1, length(n_patients)), sd = sd, mods = mods, prior_list = prior_list, n_sim = n_sim),
+  )
+})
+
+### mods param ###
+
+test_that("assessDesign validates mods parameter input and give appropriate error messages", {
+  
+  # assertions that aren't tested here for sake of brevity
+    # mods should be non-NULL object of class "Mods" from {DoseFinding}
+  
+  
+  # checking that DoseFinding didn't change how they named their 'doses' attribute
+  expect_true(
+    "doses" %in% names(attributes(mods))
+  )
+  
+  mods2 <- mods
+  attr(mods2, "doses") <- 0
+  expect_error(
+    assessDesign(n_patients = n_patients, mods = mods2, sd = sd, prior_list = prior_list, n_sim = n_sim)
+  )
+  rm(mods2)
+  
+})
+
+## prior_list param ###
+
+test_that("assessDesign validates prior_list parameter input and give appropriate error messages", {
+  
+  # assertions that aren't tested here for sake of brevity
+    # prior_list should be a non-NULL named list with length = number of dose levels
+    # length(attr(prior_list, "dose_levels")) == n_patients (see above)
+  
+  # checking that we didn't change how we named the 'dose_levels' attribute
+  expect_true(
+    "doses" %in% names(attributes(mods))
+  )
+  
+})
+
+#########################
+# Tests for getCritProb #
+#########################
+
+# getCritProb relies on DoseFinding, which we assumes works correctly, so the tests here are minimal
+
+test_that("getCritProb returns the right type of value under normal case", {
+  
+  crit_pval = getCritProb(
+    mods = mods, 
+    dose_levels = dose_levels, 
+    dose_weights = n_patients, 
+    alpha_crit_val = alpha_crit_val
+  )
+
+  expect_type(
+    crit_pval, "double"
+  )
+
+  expect_true(
+    crit_pval >= 0 & crit_pval <= 1
+  )
+
+})
+
+#########################
+# Tests for getContrMat #
+#########################
+
+# getContrMat relies on DoseFinding, which we assumes works correctly, so the tests here are minimal
+
+test_that("getContrMat returns the right type of object under normal case", {
+  
+  contr_mat = getContr(
+    mods = mods, 
+    dose_levels = dose_levels, 
+    dose_weights = n_patients,
+    prior_list = prior_list
+  )
+
+  expect_s3_class(
+    contr_mat, "optContr"
+  )
+
+})
+
+################################
+# Tests for performBayesianMCP #
+################################
+
+test_that("performBayesianMCP returns the right type of object under normal case", {
+  
+  data <- simulateData(
+    n_patients  = n_patients,
+    dose_levels = dose_levels,
+    sd          = sd,
+    mods        = mods,
+    n_sim       = n_sim
+  )
+
+  posterior_list <- getPosterior(
+    data = getModelData(data, names(mods)[1]),
+    prior_list = prior_list
+  )
+  
+  contr_mat = getContr(
+    mods = mods, 
+    dose_levels = dose_levels, 
+    dose_weights = n_patients,
+    prior_list = prior_list
+  )
+  
+  crit_pval = getCritProb(
+    mods = mods, 
+    dose_levels = dose_levels, 
+    dose_weights = n_patients, 
+    alpha_crit_val = alpha_crit_val
+  )
+
+  b_mcp <- performBayesianMCP(
+    posterior_list = posterior_list,
+    contr = contr_mat,
+    crit_prob_adj = crit_pval
+  )
+
+  expect_s3_class(
+    b_mcp,
+    "BayesianMCP"
+  )
+
+  expect_true(
+    attr(b_mcp, "crit_prob_adj") == crit_pval
+  )
+  
+  expect_type(
+    attr(b_mcp, "ess_avg"), "logical"
+  )
+
+  expect_type(
+    attr(b_mcp, "successRate"), "double"
+  )
+  
+  
+  expect_type(b_mcp, "double")
+  
+})
+
+###################################
+# Tests for performBayesianMCPMod #
+###################################
+
+test_that("performBayesianMCPMod returns the right type of object under normal case", {
+  
+  b_mcp_mod <- performBayesianMCPMod(
+    posterior_list = posterior_list,
+    contr = contr_mat,
+    crit_prob_adj = crit_pval
+  )
+
+  expect_s3_class(
+    b_mcp_mod,
+    "BayesianMCPMod"
+  )
+
+  expect_true(
+    all(names(b_mcp_mod) == c("BayesianMCP", "Mod"))
+  )
+  
+})
+
+
+#############################
+# Tests for addSignificance #
+#############################
+
+test_that("addSignificance works as intended", {
+  model_fits  <- list(linear = 1)
+
+  model_fits_with_sign = addSignificance(model_fits, list(TRUE))
+  expect_true(
+    model_fits_with_sign[[1]][["significant"]]
+  )
+
+  model_fits_with_sign = addSignificance(model_fits, list(FALSE))
+  expect_false(
+    model_fits_with_sign[[1]]$significant
+  )
+})
+
+#######################
+# Tests for BayesMCPi #
+#######################
+
+test_that("BayesMCPi function works correctly in a simple case", {
+
+  # BayesMCPi should return a list of length 3 named with "sign", "p_val", and "post_probs"
+  # The logic being tested here is: 
+    # BayesMCPi returns 1 if the posterior probability is strictly greater than the critical value, and 0 otherwise
+  
+  # Define inputs
+  posterior_i = posterior_list
+  contr_mat = list(contMat = matrix(c(0, 1), nrow = 2))
+  crit_prob = 0.5
+  
+  # Call the function
+  result = BayesMCPi(posterior_i, contr_mat, crit_prob)
+  # Check the results
+  expect_equal(result[["sign"]], 1)
+  
+  # Define inputs
+  contr_mat = list(contMat = matrix(c(0, 0), nrow = 2))
+  # Call the function
+  result = BayesMCPi(posterior_i, contr_mat, crit_prob)
+  # Check the results
+  expect_equal(result[["sign"]], rep(NA_real_, 1))
+
+#   # Define inputs
+#   contr_mat_2 = list(contMat = matrix(c(1, 0), nrow = 2))
+#   # Call the function
+#   result_2 = BayesMCPi(posterior_i, contr_mat_2, crit_prob)
+#   # Check the results
+#   expect_equal(result_2[["sign"]], 0)
+})
+
+#########################
+# Tests for getPostProb #
+#########################
+
+# Test for getPostProb
+test_that("getPostProb works correctly in a simple case", {
+  # Create a test case
+  contr_j <- c(1, 1)
+  post_combs_i <- list(
+    means = matrix(c(0, 0, 1, 1), nrow = 2),
+    vars = matrix(c(0, 0, 1, 1), nrow = 2),
+    weights = c(0.5, 0.5)
+  )
+  # Call the function with the test case
+  result <- getPostProb(
+    contr_j, 
+    post_combs_i
+  )
+  # Assert the expected behavior
+  expect_equal(result, stats::pnorm(1))
+})
diff --git a/tests/testthat/test-bootstrapping.R b/tests/testthat/test-bootstrapping.R
new file mode 100644
index 0000000..b747a0f
--- /dev/null
+++ b/tests/testthat/test-bootstrapping.R
@@ -0,0 +1,26 @@
+# Test cases
+test_that("test getBootstrapBands", {
+  
+  BMCP_result <- performBayesianMCP(
+    posterior_list = posterior_list,
+    contr           = contr_mat,
+    crit_prob_adj   = crit_pval)
+
+  model_shapes  <- c("linear", "emax", "exponential")
+  
+  # Option b) Making use of the complete posterior distribution
+  fit <- getModelFits(
+    models      = model_shapes,
+    dose_levels = dose_levels,
+    posterior   = posterior_list,
+    simple      = FALSE)
+  
+  result <- getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = c(0, 1,2,3,4,6,8))
+  expect_type(result, "list")
+  
+  result_2 <- getBootstrapQuantiles(fit, n_samples = 1e2, quantiles = c(0.1, 0.9), avg_fit = FALSE, doses = c(1, 2, 3))
+  expect_type(result_2, "list")
+
+  result_3 <- getBootstrapQuantiles(fit, quantiles = c(0.025,0.5, 0.975), doses = NULL)
+  expect_type(result_3, "list")
+})
\ No newline at end of file
diff --git a/tests/testthat/test-modelling.R b/tests/testthat/test-modelling.R
new file mode 100644
index 0000000..310668c
--- /dev/null
+++ b/tests/testthat/test-modelling.R
@@ -0,0 +1,83 @@
+# 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, se_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")
+})
+
+# Test for addModelWeights
+test_that("addModelWeights works correctly in a simple case", {
+  # Create a test case
+  test_model_fits1 <- list(
+    model1 = list(gAIC = 1),
+    model2 = list(gAIC = 1),
+    model3 = list(gAIC = 1)
+  )
+  expected_output1 <- list(
+    list(gAIC = 1, model_weight = 1/3),
+    list(gAIC = 1, model_weight = 1/3),
+    list(gAIC = 1, model_weight = 1/3)
+  )
+  
+  # Call the function with the test case
+  result <- addModelWeights(model_fits = test_model_fits1)
+  # Assert the expected behavior
+  expect_true(is.list(result))
+  expect_equal(length(result), length(test_model_fits1))
+  expect_true(all(sapply(result, function(x) "model_weight" %in% names(x))))
+  expect_equal(result, expected_output1)
+  
+  test_model_fits2 <- list(
+    model1 = list(gAIC = 100),
+    model2 = list(gAIC = 0),
+    model3 = list(gAIC = 100)
+  )
+  expected_output2 <- list(
+    list(gAIC = 100, model_weight = 0),
+    list(gAIC = 0, model_weight = 1),
+    list(gAIC = 100, model_weight = 0)
+  )
+  result <- addModelWeights(model_fits = test_model_fits2)
+  expect_equal(result, expected_output2, tolerance = 1e-10)
+})
+
+# Test for getGenAIC
+test_that("getGenAIC calculates AIC correctly using snapshot test and a simple example", {
+  # Create a test case
+  test_model_fit <- list(
+    pred_values = rep(1, 3),
+    coeffs = rep(1, 3)
+  )
+  test_post_combs <- list(
+    means = matrix(rep(1, 6), nrow = 2),
+    vars = matrix(rep(1, 6), nrow = 2),
+    weights = c(1, 1)
+  )
+  
+  result <- getGenAIC(model_fit = test_model_fit, post_combs = test_post_combs)
+  # expected result determined with snapshot of behavior prior to first release on CRAN
+  expected_result <- 6 
+  # Assert the expected behavior
+  expect_equal(result, expected_result)
+})
diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R
new file mode 100644
index 0000000..6590ffa
--- /dev/null
+++ b/tests/testthat/test-plot.R
@@ -0,0 +1,171 @@
+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(
+    posterior_list = list(posterior_emax),
+    contr          = 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
new file mode 100644
index 0000000..c26bf01
--- /dev/null
+++ b/tests/testthat/test-posterior.R
@@ -0,0 +1,107 @@
+test_that("getPosterior works correctly", {
+  dummy_data <- getModelData(data, names(mods)[1])
+  
+  # Test getPosterior function
+  posterior_list <- getPosterior(
+    data = getModelData(data, names(mods)[1]),
+    prior_list = prior_list
+  )
+  expect_type(posterior_list, "list")
+  expect_s3_class(posterior_list, "postList")
+})
+
+test_that("getPriorList input parameters do work as intented", {
+  # setup
+  library(clinDR)
+  library(dplyr)
+  set.seed(8080)
+  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)
+  
+  ##Create MAP Prior
+  hist_data <- data.frame(
+    trial = histcontrol$nctno,
+    est   = histcontrol$rslt,
+    se    = histcontrol$se,
+    sd    = histcontrol$sd,
+    n     = histcontrol$sampsize)
+  
+  dose_levels <- c(0, 2.5, 5, 10)
+  
+  # call without parameter
+  expect_error(getPriorList())
+  # only passing the historical data
+  expect_error(getPriorList(hist_data = hist_data))
+  # passing both needed parameters
+  expect_type(getPriorList(
+    hist_data = hist_data,
+    dose_levels = dose_levels,
+    robust_weight = 0.5
+  ), "list")
+  # passing wrong format for hist_data
+  expect_error(getPriorList(
+    hist_data = testdata,
+    dose_levels = dose_levels,
+    robust_weight = 0.5
+  ))
+  # passing wrong format for dose_levels
+  expect_error(getPriorList(
+    hist_data = hist_data,
+    dose_levels = c("hello", "world"),
+    robust_weight = 0.5
+  ))
+  
+})
+
+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)
+  se_hat <- matrix(c(1, 2, 3, 4), nrow = 4, ncol = 1)
+  
+  # Test getPosteriorI function
+  post_list <- getPosteriorI(data_i, prior_list, mu_hat, se_hat)
+  expect_type(post_list, "list")
+  expect_s3_class(post_list, "postList")
+  
+  # Test mu_hat and sd_hat both null branch
+  post_list <- getPosteriorI(data_i, prior_list, NULL, NULL)
+  expect_type(post_list, "list")
+  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")
+})
+
+test_that("getPostCombsI returns an object with correct attributes", {
+  posterior_i <- list(
+    matrix(c(2, 1, 2), nrow = 3),
+    matrix(c(2, 1, 2), nrow = 3)
+  )
+  result <- getPostCombsI(posterior_i)
+  
+  expect_true(is.list(result))
+  expect_equal(length(result), 3)
+  expect_equal(names(result), c("weights", "means", "vars"))
+  expect_equal(result$weights, 4)
+  expect_true(all(result$vars == c(4, 4)))
+})
diff --git a/tests/testthat/test-s3methods.R b/tests/testthat/test-s3methods.R
new file mode 100644
index 0000000..6501720
--- /dev/null
+++ b/tests/testthat/test-s3methods.R
@@ -0,0 +1,70 @@
+test_that("print.BayesianMCPMod works as intented", {
+  expect_error(print.BayesianMCPMod())
+  
+})
+
+
+test_that("print.BayesianMCP works as intented", {
+  expect_error(print.BayesianMCP())
+  
+  
+  
+})
+
+test_that("predict.ModelFits works as intented", {
+  expect_error(predict.ModelFits())
+  
+})
+
+test_that("s3 postList functions work as intented", {
+  # setup
+  library(clinDR)
+  library(dplyr)
+  set.seed(8080)
+  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)
+  
+  ##Create MAP Prior
+  hist_data <- data.frame(
+    trial = histcontrol$nctno,
+    est   = histcontrol$rslt,
+    se    = histcontrol$se,
+    sd    = histcontrol$sd,
+    n     = histcontrol$sampsize)
+  
+  dose_levels <- c(0, 2.5, 5, 10)
+  post_test_list <- getPriorList(
+    hist_data = hist_data,
+    dose_levels = dose_levels,
+    robust_weight = 0.5)
+  expect_error(summary.postList())
+  expect_type(summary.postList(post_test_list), "double")
+  expect_error(print.postList())
+  expect_type(print(post_test_list), "list")
+  expect_type(print.postList(post_test_list), "list")
+  # expect_true(names(print(post_test_list)) == c("Summary of Posterior Distributions",
+  #                                               "Maximum Difference to Control and Dose Group",
+  #                                               "Posterior Distributions"))
+})
+
+test_that("test modelFits s3 methods", {
+  
+  model_shapes <- colnames(contr_mat$contMat)
+  dose_levels  <- as.numeric(rownames(contr_mat$contMat))
+  
+  model_fits  <- getModelFits(
+    models      = model_shapes,
+    dose_levels = dose_levels,
+    posterior   = posterior_list,
+    simple      = simple)
+  
+  pred <- predict(model_fits)
+  pred_dosage <- predict(model_fits, doses = dose_levels)
+  
+  expect_type(pred, "list")
+  expect_true(is.null(attr(pred, "doses")))
+  expect_identical(attr(pred_dosage, "doses"), dose_levels)
+  expect_type(print(model_fits), "double")
+})
diff --git a/tests/testthat/test-simulation.R b/tests/testthat/test-simulation.R
new file mode 100644
index 0000000..194800c
--- /dev/null
+++ b/tests/testthat/test-simulation.R
@@ -0,0 +1,20 @@
+# Test for getModelData
+test_that("getModelData returns a data frame with correct columns", {
+
+  mock_sim_data <- data.frame(
+    simulation = 1:5,
+    dose = rnorm(5),
+    model1 = rnorm(5),
+    model2 = rnorm(5)
+  )
+  mock_model_name <- "model1"
+
+  result <- getModelData(sim_data = mock_sim_data, model_name = mock_model_name)
+  
+  # Assert that the result is a data frame
+  expect_true(is.data.frame(result))
+  # Assert that the data frame has the correct number of columns
+  expect_equal(ncol(result), 3)
+  # Assert that the data frame has the correct column names
+  expect_equal(colnames(result), c("simulation", "dose", "response"))
+})
\ No newline at end of file
diff --git a/vignettes/Simulation_Example.Rmd b/vignettes/Simulation_Example.Rmd
index 8f6c247..1ff38de 100644
--- a/vignettes/Simulation_Example.Rmd
+++ b/vignettes/Simulation_Example.Rmd
@@ -193,6 +193,4 @@ success_probabilities <- assessDesign(
   sd          = sd_tot,
   dr_means = c(-12,-14,-15,-16,-17))
 success_probabilities
-```
-
-
+```
\ No newline at end of file
diff --git a/vignettes/analysis_normal.Rmd b/vignettes/analysis_normal.Rmd
index 870c640..e3f34b3 100644
--- a/vignettes/analysis_normal.Rmd
+++ b/vignettes/analysis_normal.Rmd
@@ -207,7 +207,7 @@ contr_mat<- getContr(
 #  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}
+```{r  Execution of Bayesian MCPMod Test step}
 BMCP_result <- performBayesianMCP(
   posterior_list = posterior,
   contr          = contr_mat, 
@@ -277,6 +277,3 @@ overall_result<-performBayesianMCPMod(
 
 overall_result$Mod
 ```
-
-
-