diff --git a/NEWS.md b/NEWS.md
index fd003f0..b9218f2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,8 +1,11 @@
+# BayesSurvive (development version)
+
 # BayesSurvive 0.0.4
 
 * Add units tests
 * Rename the output of MPM coefficients in function `coef.BayesSurvive()`
 * Added `cpp` argument to `BayesSurvive()` to allow for faster computation using `Rcpp`
+* Added validation to some `BayesSurvive()` arguments
 
 # BayesSurvive 0.0.3
 
diff --git a/R/BayesSurvive.R b/R/BayesSurvive.R
index 2a9edc1..9a37fa4 100755
--- a/R/BayesSurvive.R
+++ b/R/BayesSurvive.R
@@ -111,6 +111,9 @@ BayesSurvive <- function(survObj,
                          output_graph_para = FALSE,
                          verbose = TRUE,
                          cpp = FALSE) {
+  # Validation
+  stopifnot(burnin < nIter)
+
   # same number of covariates p in all subgroups
   p <- ifelse(is.list(survObj[[1]]), NCOL(survObj[[1]]$X), NCOL(survObj$X))
   Beta.ini <- numeric(p)
diff --git a/tests/testthat/test-validation.R b/tests/testthat/test-validation.R
new file mode 100644
index 0000000..130da2f
--- /dev/null
+++ b/tests/testthat/test-validation.R
@@ -0,0 +1,25 @@
+# Load the example dataset
+dataset <- list(
+  "X" = simData[[1]]$X,
+  "t" = simData[[1]]$time,
+  "di" = simData[[1]]$status
+)
+# Run a Bayesian Cox model
+
+## Initial value: null model without covariates
+initial <- list("gamma.ini" = rep(0, ncol(dataset$X)))
+
+# Prior parameters
+hyperparPooled = list(
+  "c0"     = 2,                      # prior of baseline hazard
+  "tau"    = 0.0375,                 # sd (spike) for coefficient prior
+  "cb"     = 20,                     # sd (slab) for coefficient prior
+  "pi.ga"  = 0.02,                   # prior variable selection probability for standard Cox models
+  "a"      = -4,                     # hyperparameter in MRF prior
+  "b"      = 0.1,                    # hyperparameter in MRF prior
+  "G"      = simData$G               # hyperparameter in MRF prior
+)
+
+test_that("burnin must be less than nIter", {
+  expect_error(BayesSurvive(dataset, Iter = 30, burnin = 30))
+})