diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 906863f0..1c0ae613 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -1,6 +1,94 @@
+# ---- Library Models ----
 mod_1cmt_iv <- new_ode_model("pk_1cmt_iv")
 mod_2cmt_iv <- new_ode_model("pk_2cmt_iv")
 mod_1cmt_oral <- new_ode_model("pk_1cmt_oral")
+
+# ---- Shared ODE Models ----
+# 1-cmt oral with explicit code (used in compare tests)
+mod_1cmt_oral_code <- new_ode_model(
+  code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
+  obs = list(cmt = 2, scale = "V")
+)
+
+# Model with dose in compartment 2 (used in compare tests)
+mod_dose_cmt_2 <- new_ode_model(
+  code = "
+    dAdt[1] = -KA * A[1];
+    dAdt[2] = KA*A[1] -(CL/V) * A[2]
+    dAdt[3] = S2*(A[2]-A[3])
+  ",
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 2),
+  cpp_show_code = FALSE
+)
+
+# Basic 1-cmt without IOV (used in IOV tests)
+mod_1cmt_no_iov <- new_ode_model(
+  code = "
+    dAdt[1] = -KA * A[1]
+    dAdt[2] = +KA * A[1] -(CL/V) * A[2]
+  ",
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 1, bioav = 1),
+  parameters = c("CL", "V", "KA"),
+  cpp_show_code = FALSE
+)
+
+# 1-cmt with IOV structure (used in IOV tests)
+mod_1cmt_iov <- new_ode_model(
+  code = "
+    CL_iov = CL * exp(kappa_CL + eta_CL);
+    dAdt[1] = -KA * A[1]
+    dAdt[2] = +KA * A[1] -(CL_iov/V) * A[2]
+  ",
+  iov = list(
+    cv = list(CL = 0.2),
+    n_bins = 3
+  ),
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 1, bioav = 1),
+  declare_variables = c("kappa_CL", "CL_iov"),
+  parameters = c("kappa_CL_1", "kappa_CL_2", "kappa_CL_3", "eta_CL", "CL", "V", "KA"),
+  cpp_show_code = FALSE
+)
+
+# ---- NOT_CRAN Models (for ADVAN comparison tests) ----
+if (identical(Sys.getenv("NOT_CRAN"), "true")) {
+  # Standard parameters for ADVAN tests
+  advan_parameters <- list(
+    CL = 10, V = 50, KA = 0.5,
+    Q = 5, V2 = 100,
+    Q2 = 3, V3 = 150,
+    F1 = 1
+  )
+
+  # 1-cmt IV with AUC compartment
+  mod_1cmt_iv_auc <- new_ode_model(
+    code = "dAdt[1] = -(CL/V)*A[1]; dAdt[2] = A[1]/V;",
+    parameters = advan_parameters
+  )
+
+  # 2-cmt IV with AUC compartment
+  mod_2cmt_iv_auc <- new_ode_model(
+    code = "
+      dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2];
+      dAdt[2] = +(Q/V)*A[1] - (Q/V2)*A[2];
+      dAdt[3] = A[1]/V;
+    ",
+    parameters = advan_parameters
+  )
+
+  # 3-cmt IV with AUC compartment
+  mod_3cmt_iv_auc <- new_ode_model(
+    code = "
+      dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2] - (Q2/V)*A[1] + (Q2/V3)*A[3];
+      dAdt[2] =                (Q/V)*A[1]  -(Q/V2)*A[2]                             ;
+      dAdt[3] =                                           (Q2/V)*A[1] - (Q2/V3)*A[3];
+      dAdt[4] = A[1]/V;
+    ",
+    parameters = advan_parameters
+  )
+}
 mod_1cmt_oral_lagtime <- new_ode_model(
   code = "
     dAdt[0] = -KA * A[0]
diff --git a/tests/testthat/test-advan.R b/tests/testthat/test-advan.R
new file mode 100644
index 00000000..a2e95cb9
--- /dev/null
+++ b/tests/testthat/test-advan.R
@@ -0,0 +1,325 @@
+# ---- Shared Setup ----
+## These models are also tested in the unit tests for `calc_ss_analytics()`, so just testing a few example cases here
+dose <- 100
+interval <- 12
+t_inf <- 1
+n_days <- 5
+parameters <- list(CL = 10, V = 50, KA = 0.5, Q = 5, V2 = 100, Q2 = 3, V3 = 150, F1 = 1)
+t_obs <- c(3, 6, 8, 23)
+reg_bolus <- new_regimen(
+  amt = dose,
+  times = seq(0, interval * n_days * (24/interval), interval),
+  t_inf = t_inf, type = "bolus"
+)
+data <- advan_create_data(
+  reg_bolus,
+  parameters = parameters,
+  cmts = 5,
+  t_obs = t_obs
+)
+
+## Infusion dataset
+reg_infusion <- new_regimen(
+  amt = dose,
+  times = seq(0, interval * n_days * (24/interval), interval),
+  t_inf = t_inf,
+  type = "infusion"
+)
+data_infusion <- advan_create_data(
+  reg_infusion,
+  parameters = parameters,
+  cmts = 6,
+  t_obs = t_obs
+)
+
+# ---- One Compartment Tests ----
+test_that("One compartment IV bolus", {
+  res1_iv   <- advan("1cmt_iv_bolus", cpp=FALSE)(data)
+  res1_iv_c <- advan("1cmt_iv_bolus", cpp=TRUE)(data)
+  expect_equal(round(res1_iv[res1_iv$TIME == 23,]$DV, 3), 0.242)
+  expect_true(!any(is.na(res1_iv$DV)))
+  expect_equal(res1_iv, res1_iv_c)
+})
+
+test_that("One compartment IV infusion", {
+  res1_iv_inf   <- advan("1cmt_iv_infusion", cpp=FALSE)(data_infusion)
+  res1_iv_inf_c <- advan("1cmt_iv_infusion", cpp=TRUE)(data_infusion)
+  f1 <- advan("1cmt_iv_infusion", cpp=FALSE)
+  f2 <- advan("1cmt_iv_infusion", cpp=TRUE)
+
+  expect_equal(round(res1_iv_inf[res1_iv_inf$TIME == 23,]$DV, 3), 0.268)
+  expect_true(!any(is.na(res1_iv_inf$DV)))
+  expect_equal(res1_iv_inf, res1_iv_inf_c)
+  expect_equal(attr(f1, "type"), "infusion")
+  expect_equal(attr(f2, "type"), "infusion")
+  expect_equal(attr(f1, "implementation"), FALSE)
+  expect_equal(attr(f2, "implementation"), TRUE)
+  expect_equal(attr(f1, "cmt"), 1)
+  expect_equal(attr(f2, "cmt"), 1)
+})
+
+test_that("One compartment oral", {
+  res1_oral <- advan("1cmt_oral", cpp=FALSE)(data)
+  res1_oral_c <- advan("1cmt_oral", cpp=TRUE)(data)
+
+  expect_equal(round(res1_oral[res1_oral$TIME == 23,]$DV, 3), 0.389)
+  expect_true(!any(is.na(res1_oral$DV)))
+  expect_equal(res1_oral, res1_oral_c)
+})
+
+# ---- Two Compartment Tests ----
+test_that("Two compartment iv bolus", {
+  res2_iv   <- advan("2cmt_iv_bolus", cpp=FALSE)(data)
+  res2_iv_c <- advan("2cmt_iv_bolus", cpp=TRUE)(data)
+
+  expect_equal(round(res2_iv[res2_iv$TIME == 23,]$DV, 3), 0.212)
+  expect_true(!any(is.na(res2_iv$DV)))
+  expect_equal(res2_iv, res2_iv_c)
+})
+
+test_that("Two compartment iv infusion", {
+  res2_iv_inf   <- advan("2cmt_iv_infusion", cpp=FALSE)(data_infusion)
+  res2_iv_inf_c <- advan("2cmt_iv_infusion", cpp=TRUE)(data_infusion)
+
+  expect_equal(round(res2_iv_inf[res2_iv_inf$TIME == 23,]$DV, 3), 0.225)
+  expect_true(!any(is.na(res2_iv_inf$DV)))
+  expect_equal(res2_iv_inf, res2_iv_inf_c)
+})
+
+test_that("Two compartment oral", {
+  res2_oral   <- advan("2cmt_oral", cpp=FALSE)(data)
+  res2_oral_c <- advan("2cmt_oral", cpp=TRUE)(data)
+
+  expect_equal(round(res2_oral[res2_oral$TIME == 23,]$DV, 3), 0.302)
+  expect_true(!any(is.na(res2_oral$DV)))
+  expect_equal(res2_oral, res2_oral_c)
+})
+
+# ---- Three Compartment Tests ----
+test_that("Three compartment IV bolus", {
+  res3_iv   <- advan("3cmt_iv_bolus", cpp=FALSE)(data)
+  res3_iv_c <- advan("3cmt_iv_bolus", cpp=TRUE)(data)
+
+  expect_equal(round(res3_iv[res3_iv$TIME == 23,]$DV, 3), 0.169)
+  expect_true(!any(is.na(res3_iv$DV)))
+  expect_equal(res3_iv, res3_iv_c)
+})
+
+test_that("Three compartment IV infusion", {
+  res3_iv_inf   <- advan("3cmt_iv_infusion", cpp=FALSE)(data_infusion)
+  res3_iv_inf_c <- advan("3cmt_iv_infusion", cpp=TRUE)(data_infusion)
+
+  expect_equal(round(res3_iv_inf[res3_iv_inf$TIME == 23,]$DV, 3), 0.177)
+  expect_true(!any(is.na(res3_iv_inf$DV)))
+  expect_equal(res3_iv_inf, res3_iv_inf_c)
+})
+
+test_that("Three compartment IV oral", {
+  res3_oral <- advan("3cmt_oral", cpp=FALSE)(data)
+  res3_oral_c <- advan("3cmt_oral", cpp=TRUE)(data)
+
+  expect_equal(round(res3_oral[res3_oral$TIME == 23,]$DV, 3), 0.236)
+  expect_true(!any(is.na(res3_oral$DV)))
+  expect_equal(res3_oral, res3_oral_c)
+})
+
+# ---- AUC Comparison Tests (NOT_CRAN) ----
+if (identical(Sys.getenv("NOT_CRAN"), "true")) {
+  # Setup for AUC tests
+  t_obs_auc <- c(3, 6, 8, 23, 47)
+  t_inf_auc <- 1.5
+
+  ## bolus dataset for AUC tests
+  reg_bolus_auc <- new_regimen(
+    amt = dose,
+    times = seq(0, interval * n_days * (24/interval), interval),
+    t_inf = t_inf_auc,
+    type = "bolus"
+  )
+  data_bolus_auc <- advan_create_data(
+    reg_bolus_auc,
+    parameters = parameters,
+    cmts = 5,
+    t_obs = t_obs_auc
+  )
+
+  ## Infusion dataset for AUC tests
+  reg_infusion_auc <- new_regimen(
+    amt = dose,
+    times = seq(0, interval * n_days * (24/interval), interval),
+    t_inf = t_inf_auc,
+    type = "infusion"
+  )
+  data_infusion_auc <- advan_create_data(
+    reg_infusion_auc,
+    parameters = parameters,
+    cmts = 6,
+    t_obs = t_obs_auc
+  )
+}
+
+test_that("One compartment bolus ADVAN runs", {
+  skip_on_cran()
+  res1_iv_r <- advan("1cmt_iv_bolus", cpp=FALSE)(data_bolus_auc)
+  res1_iv_c <- advan("1cmt_iv_bolus", cpp=TRUE)(data_bolus_auc)
+  res1_iv_ode <- sim(ode = mod_1cmt_iv_auc, regimen = reg_bolus_auc, parameters = parameters, t_obs = t_obs_auc)
+
+  # AUC R
+  expect_equal(round(res1_iv_r[res1_iv_r$TIME %in% t_obs_auc,]$AUC, 5), round(res1_iv_ode[res1_iv_ode$comp == 2,]$y, 5))
+
+  #AUC-C
+  expect_equal(round(res1_iv_c[res1_iv_c$TIME %in% t_obs_auc,]$AUC, 5), round(res1_iv_ode[res1_iv_ode$comp == 2,]$y, 5))
+})
+
+test_that("Two compartment bolus ADVAN runs", {
+  skip_on_cran()
+  res2_iv_r <- advan("2cmt_iv_bolus", cpp=FALSE)(data_bolus_auc)
+  res2_iv_c <- advan("2cmt_iv_bolus", cpp=TRUE)(data_bolus_auc)
+  res2_iv_ode <- sim(ode = mod_2cmt_iv_auc, regimen = reg_bolus_auc, parameters = parameters, t_obs = t_obs_auc)
+  expect_equal(
+    round(res2_iv_r[res2_iv_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res2_iv_c[res2_iv_c$TIME %in% t_obs_auc,]$AUC, 5)
+  )
+  # AUC R
+  expect_equal(
+    round(res2_iv_r[res2_iv_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res2_iv_ode[res2_iv_ode$comp == 3,]$y, 5)
+  )
+
+  #AUC-C
+  expect_equal(
+    round(res2_iv_c[res2_iv_c$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res2_iv_ode[res2_iv_ode$comp == 3,]$y, 5)
+  )
+})
+
+test_that("Two compartment infusion ADVAN runs", {
+  skip_on_cran()
+  res2_inf_r <- advan("2cmt_iv_infusion", cpp=FALSE)(data_infusion_auc)
+  res2_inf_c <- advan("2cmt_iv_infusion", cpp=TRUE)(data_infusion_auc)
+  res2_inf_ode <- sim(ode = mod_2cmt_iv_auc, regimen = reg_infusion_auc, parameters = parameters, t_obs = t_obs_auc)
+
+  expect_equal(
+    round(res2_inf_r[res2_inf_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res2_inf_c[res2_inf_c$TIME %in% t_obs_auc,]$AUC, 5)
+  )
+
+  # AUC R
+  expect_equal(
+    round(res2_inf_r[res2_inf_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res2_inf_ode[res2_inf_ode$comp == 3,]$y, 5)
+  )
+
+  #AUC-C
+  expect_equal(
+    round(res2_inf_c[res2_inf_c$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res2_inf_ode[res2_inf_ode$comp == 3,]$y, 5)
+  )
+
+})
+
+test_that("Three compartment bolus ADVAN runs", {
+  skip_on_cran()
+  res3_iv_r <- advan("3cmt_iv_bolus", cpp=FALSE)(data_bolus_auc)
+  res3_iv_c <- advan("3cmt_iv_bolus", cpp=TRUE)(data_bolus_auc)
+  res3_iv_ode <- sim(ode = mod_3cmt_iv_auc, regimen = reg_bolus_auc, parameters = parameters, t_obs = t_obs_auc)
+  expect_equal(
+    round(res3_iv_r[res3_iv_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res3_iv_c[res3_iv_c$TIME %in% t_obs_auc,]$AUC, 5)
+  )
+  # AUC R
+  expect_equal(
+    round(res3_iv_r[res3_iv_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
+  )
+
+  #AUC-C
+  expect_equal(
+    round(res3_iv_c[res3_iv_c$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
+  )
+})
+
+test_that("Three compartment iv ADVAN runs", {
+  skip_on_cran()
+  res3_iv_r <- advan("3cmt_iv_infusion", cpp=FALSE)(data_infusion_auc)
+  res3_iv_c <- advan("3cmt_iv_infusion", cpp=TRUE)(data_infusion_auc)
+  res3_iv_ode <- sim(ode = mod_3cmt_iv_auc, regimen = reg_infusion_auc, parameters = parameters, t_obs = t_obs_auc)
+  expect_equal(
+    round(res3_iv_r[res3_iv_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res3_iv_c[res3_iv_c$TIME %in% t_obs_auc,]$AUC, 5)
+  )
+  # AUC R
+  expect_equal(
+    round(res3_iv_r[res3_iv_r$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
+  )
+
+  #AUC-C
+  expect_equal(
+    round(res3_iv_c[res3_iv_c$TIME %in% t_obs_auc,]$AUC, 5),
+    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
+  )
+})
+
+# ---- Covariate Tests ----
+test_that("Analytic and ODE models with covariates are the same", {
+  skip_on_cran()
+
+  ## Create dataset
+  dose <- 100
+  interval <- 12
+  n_days <- 2
+  parameters <- list(
+    CL = 10,
+    V = 50,
+    KA = 0.5,
+    Q = 5,
+    V2 = 100,
+    Q2 = 3,
+    V3 = 150,
+    F1 = 1
+  )
+  t_obs <- seq(0, 40, .1)
+  reg_bolus <- new_regimen(
+    amt = dose,
+    times = seq(0, interval * n_days * (24/interval), interval),
+    type = "bolus"
+  )
+  ## there is slight difference in how bolus doses are handled.
+  ## Analytical equation is perhaps more consistent, so not testing
+  ## simulations at dose times. Should look into later.
+  t_obs <- t_obs[! t_obs %in% reg_bolus$dose_times]
+  covariates <- list(WT = new_covariate(80), CRCL=new_covariate(4.5))
+
+  ## Using analytic equations model:
+  data_ana <- sim(
+    analytical = "1cmt_iv_bolus",
+    parameters = parameters,
+    covariates = covariates,
+    regimen = reg_bolus,
+    t_obs = t_obs,
+    covariate_model = "CL = CL * (CRCL / 3)^0.75; V = V * (WT / 70.0)"
+  )
+
+  ## Using ODE model:
+  mod1 <- new_ode_model(
+    code = "
+    dAdt[1] = -( (CL*pow(CRCL/3.0, 0.75)) / (V*WT/70.0) ) * A[1];
+  ",
+  covariates = covariates,
+  obs = list(cmt = 1, scale = "V*WT/70.0"), dose = list(cmt = 1)
+  )
+  data_ode <- sim(
+    ode = mod1,
+    parameters = parameters,
+    covariates = covariates,
+    regimen = reg_bolus,
+    t_obs = t_obs,
+    duplicate_t_obs = TRUE,
+    only_obs = TRUE
+  )
+
+  expect_equal(nrow(data_ana), nrow(data_ode))
+  expect_equal(round(data_ana$y,4), round(data_ode$y, 4))
+})
diff --git a/tests/testthat/test-regimen.R b/tests/testthat/test-regimen.R
new file mode 100644
index 00000000..fc4631ee
--- /dev/null
+++ b/tests/testthat/test-regimen.R
@@ -0,0 +1,285 @@
+# ---- new_regimen Tests ----
+test_that("Basic regimen creation working", {
+  reg <- new_regimen(amt=100, n=4, interval=4)
+  expect_true(all(c("regimen", "list") %in% class(reg)))
+})
+
+test_that("Regimen type parsed correctly", {
+  reg_oral <- new_regimen(amt = 100, n = 4, interval = 4, type = "oral")
+  reg_bolus <- new_regimen(amt = 100, n=4, interval = 4, type = "bolus")
+  reg_inf <- new_regimen(amt = 100, n = 4, interval = 4, type = "infusion")
+  reg_sc <- new_regimen(amt = 100, n = 4, interval = 4, type = "sc")
+  reg_im <- new_regimen(amt = 100, n = 4, interval = 4, type = "im")
+  reg_mixed <- new_regimen(amt = 100, n = 5, interval = 4, type = c("bolus", "infusion", "oral", "sc", "im"))
+  reg_oral_susp <- new_regimen(amt = 1, n = 4, interval = 8, type = "oral_susp")
+
+  expect_equal(reg_oral$type, rep("oral", 4))
+  expect_equal(reg_bolus$type, rep("bolus", 4))
+  expect_equal(reg_inf$type, rep("infusion", 4))
+  expect_equal(reg_sc$type, rep("sc", 4))
+  expect_equal(reg_im$type, rep("im", 4))
+  expect_equal(reg_mixed$type, c("bolus", "infusion", "oral", "sc", "im"))
+  expect_equal(reg_oral_susp$type, rep("oral_susp", 4))
+
+  # oral-type regimens should have t_inf set to 0
+  expect_equal(reg_oral$t_inf, rep(0, 4))
+  expect_equal(reg_oral_susp$t_inf, rep(0, 4))
+})
+
+test_that("Auto-detect infusion vs bolus", {
+  expect_silent(reg1 <- new_regimen(amt = 100, n = 4, interval = 4, t_inf = 0))
+  expect_warning(reg2 <- new_regimen(amt = 100, n = 4, interval = 4, t_inf = 1))
+  expect_warning(reg3 <- new_regimen(amt = 100, n = 4, interval = 4, t_inf = c(0, 1, 0, 1)))
+
+  expect_equal(reg1$type, rep("bolus", 4))
+  expect_equal(reg2$type, rep("infusion", 4))
+  expect_equal(reg3$type, c("bolus", "infusion", "bolus", "infusion"))
+})
+
+test_that("n = 0 doses does not create doses at negative times", {
+  reg0 <- new_regimen(amt = 200, n = 0, interval = 24, t_inf = 2, type = "infusion")
+  expect_false(min(reg0$dose_times) < 0)
+})
+
+test_that("Rate argument creates valid PKPD regimen", {
+  reg1 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36, 48),
+    type = "infusion",
+    rate = c(1,2,3,4,5)
+  )
+  expect_equal(round(reg1$t_inf), c(100, 50, 33, 25, 20))
+
+  reg2 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36, 48),
+    type = "infusion",
+    rate = c(5)
+  )
+  expect_equal(round(reg2$t_inf), rep(20, 5))
+})
+
+test_that("Doses < 0 set to 0", {
+  expect_warning(
+    tmp <- new_regimen(amt = c(-1, -2, 3, 4), times = c(0, 24, 48, 72), type = "infusion")
+  )
+  expect_true(all(tmp$dose_amts >= 0))
+})
+
+test_that("new_regimen can take arbitrary values for `type`", {
+  reg <- new_regimen(100, times = 0, type = "pip")
+  expect_equal(reg$type, "pip")
+})
+
+test_that("do not creat regimens of `type` 'covariate'", {
+  expect_error(new_regimen(100, times = 0, type = "covariate"))
+})
+
+test_that("sc doses accept an infusion length argument'", {
+  reg1 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36, 48),
+    type = "sc",
+    t_inf = 30/60
+  )
+  expect_equal(reg1$t_inf, rep(0.5,5))
+})
+
+test_that("t_inf imputed correctly", {
+  reg1 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36, 48, 60, 72, 84),
+    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral")
+  )
+  reg2 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36, 48, 60, 72, 84),
+    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral"),
+    t_inf = c(2/60, 2.5, 3/60, NA, NA, NA, NA, NA)
+  )
+  reg3 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36, 48, 60, 72, 84),
+    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral"),
+    t_inf = numeric(0)
+  )
+  reg4 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36, 48, 60, 72, 84),
+    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral"),
+    t_inf = NULL
+  )
+  reg5 <- new_regimen(
+    amt = 100,
+    times = c(0, 12, 24, 36),
+    type = c("sc", "infusion", "im", "unknown_drug_type"),
+    t_inf = c(2/60, 2.5, 3/60, NA)
+  )
+  expect_equal(reg1$t_inf, c(1/60, 1, 1/60, 1/60, 1, 1/60, 0, 0))
+  expect_equal(reg2$t_inf, c(2/60, 2.5, 3/60, 1/60, 1, 1/60, 0, 0))
+  expect_equal(reg3$t_inf, c(1/60, 1, 1/60, 1/60, 1, 1/60, 0, 0))
+  expect_equal(reg4$t_inf, c(1/60, 1, 1/60, 1/60, 1, 1/60, 0, 0))
+  expect_equal(reg5$t_inf, c(2/60, 2.5, 3/60, 1))
+})
+
+# ---- merge_regimen Tests ----
+test_that("merge regimens correctly", {
+  reg1 <- new_regimen(amt = 1000, times = c(0, 24), type = "infusion", t_inf = 1)
+  reg2 <- new_regimen(amt = 500, times = c(12, 36), type = "oral")
+  reg3 <- merge_regimen(regimens = list(reg1, reg2))
+  expect_equal(reg3$dose_times, c(0, 12, 24, 36))
+  expect_equal(reg3$type, c("infusion", "oral", "infusion", "oral"))
+  expect_equal(reg3$dose_amts, c(1000, 500, 1000, 500))
+  expect_equal(reg3$t_inf, c(1, 0, 1, 0))
+  expect_null(reg3$cmt)
+})
+
+test_that("merge regimens correctly: two infusions", {
+  reg1 <- new_regimen(amt = 1000, times = c(0, 24), type = "infusion_1", t_inf = 1)
+  reg2 <- new_regimen(amt = 500, times = c(12, 36), type = "infusion_2", t_inf = 2)
+  reg3 <- merge_regimen(regimens = list(reg1, reg2))
+  expect_equal(reg3$dose_times, c(0, 12, 24, 36))
+  expect_equal(reg3$type, c("infusion_1", "infusion_2", "infusion_1", "infusion_2"))
+  expect_equal(reg3$dose_amts, c(1000, 500, 1000, 500))
+  expect_equal(reg3$t_inf, c(1, 2, 1, 2))
+  expect_null(reg3$cmt)
+})
+
+test_that("cmt info gets merged as well, when available", {
+  reg1 <- new_regimen(amt = 1000, times = c(0, 24), type = "infusion", t_inf = 1, cmt = 2)
+  reg2 <- new_regimen(amt = 500, times = c(12, 36), type = "oral", cmt = 1)
+  reg3 <- merge_regimen(regimens = list(reg1, reg2))
+  expect_equal(reg3$dose_times, c(0, 12, 24, 36))
+  expect_equal(reg3$type, c("infusion", "oral", "infusion", "oral"))
+  expect_equal(reg3$dose_amts, c(1000, 500, 1000, 500))
+  expect_equal(reg3$t_inf, c(1, 0, 1, 0))
+  expect_equal(reg3$cmt, c(2, 1, 2, 1))
+})
+
+# ---- join_regimen Tests ----
+# Shared regimens for join tests
+reg1_join <- new_regimen(
+  amt = c(100, 100, 100),
+  t_inf = c(24, 24, 28),
+  time = c(0, 24, 48),
+  type = "infusion"
+)
+reg2_join <- new_regimen(
+  amt = c(200, 200, 200),
+  t_inf = c(24,24,24),
+  time = c(0, 24, 48),
+  type = "infusion"
+)
+
+test_that("Early update of a regimen works", {
+  res1 <- join_regimen(reg1_join, reg2_join, t_dose_update = 65)
+  res2 <- join_regimen(reg1_join, reg2_join, interval = 15)
+  expect_equal(res1$t_inf, c(24, 24, 17, 24, 24, 24))
+  expect_equal(round(res1$dose_amts[3], 1), round(100 * 17 / 28, 1))
+  expect_equal(res2$t_inf, c(24, 24, 15, 24, 24, 24))
+  expect_equal(round(res2$dose_amts[3], 1), round(100 * 15 / 28, 1))
+})
+
+test_that("Gap in regimen is ok", {
+  res3 <- join_regimen(reg1_join, reg2_join, t_dose_update = 100)
+  expect_equal(res3$t_inf, c(24, 24, 28, 24, 24, 24))
+  expect_equal(res3$dose_amts, c(100, 100, 100, 200, 200, 200))
+})
+
+test_that("Update from time = 0", {
+  res4 <- join_regimen(reg1_join, reg2_join, t_dose_update = 0)
+  expect_equal(res4$dose_amts, reg2_join$dose_amts)
+})
+
+test_that("Early update of a regimen works when mixed routes", {
+  reg3 <- new_regimen(
+    amt = c(100, 100, 100),
+    t_inf = c(0, 0, 0),
+    time = c(0, 24, 48),
+    type = "oral"
+  )
+  res0 <- join_regimen(reg3, reg2_join, dose_update = 1)
+  expect_equal(res0$type, rep("infusion", 3))
+
+  res1 <- join_regimen(reg3, reg2_join, dose_update = 3)
+  expect_equal(res1$type, c(rep("oral", 2), rep("infusion", 3)))
+
+  res2 <- join_regimen(reg3, reg2_join, dose_update = 4, interval = 24)
+  expect_equal(res2$type, c(rep("oral", 3), rep("infusion", 3)))
+})
+
+test_that("insufficient join timing throws errors", {
+  expect_error(
+    join_regimen(reg1_join, reg2_join)
+  )
+})
+
+test_that("dose_update longer than reg 1 length without interval specified throws error", {
+  expect_error(
+    join_regimen(reg1_join, reg2_join, dose_update = 4)
+  )
+})
+
+test_that("one null regimen returns the other", {
+  expect_equal(join_regimen(NULL, reg2_join), reg2_join)
+  expect_equal(join_regimen(reg1_join, NULL), reg1_join)
+})
+
+test_that("join regimen works when no t_inf specified, e.g. oral or sc dosing", {
+  reg_i <- new_regimen(amt = 500, type = "sc", times = c(0, 2, 6)*24*7)
+  reg_m <- new_regimen(amt = 400, type = "sc", n = 5, interval = 8*7*24)
+  reg <- join_regimen(reg_i, reg_m, interval = 8*7*24)
+  expect_equal(reg$dose_times, c(0, 336, 1008, 2352, 3696, 5040, 6384, 7728))
+  expect_equal(reg$type, c("sc", "sc", "sc", "sc", "sc", "sc", "sc", "sc"))
+  expect_equal(reg$dose_amts, c(500, 500, 500, 400, 400, 400, 400, 400))
+})
+
+test_that("join_regimen with t_dose_update maintains correct rate length", {
+  reg1 <- new_regimen(
+    amt = 1,
+    times = c(0, 24),
+    type = "oral"
+  )
+  reg2 <- new_regimen(
+    amt = 2,
+    time = c(0, 24, 48),
+    type = "oral"
+  )
+  reg3 <- join_regimen(reg1, reg2, t_dose_update = 24)
+
+  expect_equal(length(reg3$rate), length(reg3$dose_times))
+  expect_equal(reg3$n, 4)
+})
+
+# ---- shift_regimen Tests ----
+test_that("Shifting by 1 works", {
+  keys <- c("oral", "infusion", "bolus")
+  for(key in keys) {
+    reg1 <- new_regimen(amt = 2000, interval = 24, n = 6, type = key, t_inf = 1)
+    reg1_s1 <- shift_regimen(reg1)
+    expect_true("regimen" %in% class(reg1_s1))
+    expect_equal(length(reg1_s1$dose_amts), 5)
+  }
+})
+
+test_that("Shifting by N works", {
+  keys <- c("oral", "infusion", "bolus")
+  for(key in keys) {
+    reg2 <- new_regimen(amt = c(1:6), times = c(0:5) * 24, type = "infusion")
+    reg2_s1 <- shift_regimen(reg2, n = 3)
+    expect_true("regimen" %in% class(reg2_s1))
+    expect_equal(length(reg2_s1$dose_amts), 3)
+    expect_equal(length(reg2_s1$t_inf), 3)
+    # amounts taken from end, not front
+    expect_equal(reg2_s1$dose_amts, c(4, 5, 6))
+  }
+})
+
+test_that("Shifting by N > length(regimen) produces NULL", {
+  keys <- c("oral", "infusion", "bolus")
+  for(key in keys) {
+    reg2 <- new_regimen(amt = c(1:6), times = c(0:5) * 24, type = "infusion")
+    reg2_s2 <- shift_regimen(reg2, n = 10)
+    expect_null(reg2_s2)
+  }
+})
diff --git a/tests/testthat/test_sim.R b/tests/testthat/test-sim.R
similarity index 64%
rename from tests/testthat/test_sim.R
rename to tests/testthat/test-sim.R
index cb676754..ad6725d4 100644
--- a/tests/testthat/test_sim.R
+++ b/tests/testthat/test-sim.R
@@ -1,3 +1,4 @@
+# ---- Shared Setup ----
 reg <- new_regimen(
   amt = 100,
   n = 3,
@@ -10,6 +11,7 @@ par <- list(CL = 5, V = 50)
 omega <- c(0.1, 0.0, 0.1)
 t_obs <- c(2, 48)
 
+# ---- Event Table Tests ----
 test_that("return_event_table=TRUE returns an appropriate event table", {
   evtab1 <- sim_ode(
     mod_1cmt_iv,
@@ -72,6 +74,7 @@ test_that("return_event_table=TRUE returns an appropriate event table with covar
   ))
 })
 
+# ---- Bioavailability Tests ----
 test_that("sim works properly for a model where bioavailability is dependent on dose", {
   skip_on_cran()
   reg <- new_regimen(amt = 1000, n = 4, interval = 12, type = 'oral')
@@ -116,6 +119,7 @@ test_that("sim works properly for a model where bioavailability is dependent on
   expect_true(all(round(dat$F1_avg, 5) == 1)) # should all be 1 from the first value, with no NAs
 })
 
+# ---- Covariate and t_init Tests ----
 test_that("covariates and doses are shifted correctly when t_init != 0", {
   # example covs use both interpolation and locf methods.
   covs <- list(
@@ -219,7 +223,7 @@ test_that("covariates_table and doses are shifted correctly when t_init != 0", {
   expect_equal(first_dose, evtab1[which(evtab1$dose > 0)[1], ])
 })
 
-
+# ---- Steady State Covariate Tests ----
 test_that("covariates are shifted correctly when t_ss != 0", {
   covs <- list(
     SCR = new_covariate(
@@ -531,6 +535,7 @@ test_that("times are recalculated correctly after steady-state regimen added", {
   )
 })
 
+# ---- t_max Tests ----
 test_that("t_max is shifted correctly when t_ss != 0", {
   reg <- new_regimen(
     amt = 1000, interval = 12, n = 6, t_inf = 1, type = "infusion",
@@ -578,3 +583,329 @@ test_that("t_max is shifted correctly when t_ss != 0", {
   ))
   expect_equal(sum(is.na(evtab3)), 0)
 })
+
+# ---- sim_core Tests ----
+test_that("sim core works", {
+  reg <- new_regimen(amt = 100, n = 5, interval = 12, t_inf = 1, type = "infusion")
+  par <- list(KA = 1, CL = 5, V = 50)
+
+  ## have to be explicit about t_obs with sim_core!
+  f1 <- function() {
+    res <- sim(ode = mod_1cmt_oral, regimen = reg, parameters = par, only_obs = TRUE, t_obs=c(0:24))$y
+    return(res)
+  }
+
+  f2 <- function() {
+    obj <- sim(ode = mod_1cmt_oral, regimen = reg, parameters = par, only_obs = TRUE, t_obs=c(0:24), return_design=TRUE)
+    sim_core(obj, ode = mod_1cmt_oral)$y
+  }
+  expect_equal(f1(), f2())
+})
+
+test_that("sim core works for absorption model with lagtime", {
+  reg <- new_regimen(amt = 100, n = 5, interval = 12, t_inf = 1, type = "infusion")
+  par <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
+
+  ## have to be explicit about t_obs with sim_core!
+  f_ref <- function() {
+    res <- sim(
+      ode = mod_1cmt_oral_lagtime,
+      regimen = reg,
+      parameters = par,
+      only_obs = TRUE,
+      t_obs=c(0:24)
+    )$y
+    return(res)
+  }
+  f_core <- function() {
+    obj <- sim(
+      ode = mod_1cmt_oral_lagtime,
+      regimen = reg,
+      parameters = par,
+      only_obs = TRUE,
+      t_obs=c(0:24),
+      return_design=TRUE
+    )
+    sim_core(
+      obj,
+      ode = mod_1cmt_oral_lagtime,
+      lagtime = c(0.83, 0)
+    )$y ## Lagtime parameter needed!
+  }
+  expect_equal(f_ref(), f_core())
+})
+
+# ---- Lagtime Tests ----
+test_that("dose dump after lagtime in correct order in output data", {
+  skip_on_cran()
+  reg <- new_regimen(amt = 500, n = 4, interval = 12, type = 'oral')
+  pars <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
+  dat <- sim_ode(
+    ode = mod_1cmt_oral_lagtime,
+    regimen = reg,
+    parameters = pars,
+    only_obs = FALSE
+  )
+  ## Change after RXR-2394: time of TLAG not in dataset anymore unless requested by user in t_obs
+  ## Before change: expect_equal(round(dat[dat$t == 12.83 & dat$comp == 1,]$y, 1), c(1.2, 501.2))
+  ## After change:
+  expect_equal(nrow(dat[dat$t == 12.83,]), 0)
+  ## When grid requested by user, lagtime should be visible
+  dat <- sim_ode(
+    ode = mod_1cmt_oral_lagtime,
+    regimen = reg,
+    parameters = pars,
+    t_obs = seq(0, 1, 0.01),
+    only_obs = TRUE
+  )
+  tmp <- dat[dat$t >= 0.82 & dat$t <= 0.85, ]
+  expect_equal(tmp$t, c(0.82, 0.83, 0.84, 0.85))
+  expect_equal(round(tmp$y, 2), c(0, 0, 0.05, 0.10))
+})
+
+# ---- Result Comparison Tests ----
+# Analytic solution for 1-cmt oral
+pk1cmt_oral_anal <- function(t, dose, KA, V, CL) {
+  dose*KA/(V*(KA-CL/V))*(exp(-(CL/V) * t)-exp(-KA * t))
+}
+
+test_that("Library and custom C++ and code matches analytic soln", {
+  p <- list(KA = 1, CL = 5, V = 50)
+  t_obs <- c(0:72)
+  t_obs2 <- t_obs + 0.1234 # also needs to be producing results with non-integer times
+  dose <- 100
+  t_dose <- c(0)
+  regimen <- new_regimen(amt=dose, times = t_dose, type = "oral")
+
+  pk1cmt_oral_lib <- sim_ode(
+    ode = mod_1cmt_oral,
+    parameters = p,
+    regimen = regimen,
+    t_obs = t_obs,
+    int_step_size = 0.1,
+    duplicate_t_obs = TRUE,
+    only_obs=TRUE
+  )
+
+  pk1cmt_oral_code_res <- sim_ode(
+    ode = mod_1cmt_oral_code,
+    parameters = p,
+    duplicate_t_obs = TRUE,
+    regimen=regimen,
+    t_obs=t_obs,
+    int_step_size = 0.1,
+    only_obs=TRUE
+  )
+
+  pk1cmt_oral_anal_res <- pk1cmt_oral_anal(t_obs, dose, p$KA, p$V, p$CL)
+  expect_equal(round(pk1cmt_oral_lib$y, 3), round(pk1cmt_oral_anal_res, 3))
+  expect_equal(round(pk1cmt_oral_code_res$y, 3), round(pk1cmt_oral_anal_res, 3))
+})
+
+test_that("precision in time does not impact # obs returned", {
+  regimen_mult <- new_regimen(
+    amt = rep(12.8, 3),
+    times = c(0, 6, 12),
+    type="infusion",
+    t_inf = 2
+  )
+  t_obs <- c(11.916, 14.000, 16.000, 17.000, 30)
+  tmp <- sim_ode(
+    ode = mod_1cmt_iv,
+    parameters = list(CL = 5, V = 50),
+    regimen = regimen_mult,
+    t_obs = t_obs,
+    only_obs = TRUE
+  )
+  expect_equal(tmp$t, t_obs)
+})
+
+test_that("test bug EmCo 20150925", {
+  xtim <- c(0, 2, 4, 8, 12, 24)
+  sujdos <- 320
+  param <- list(KA = 1.8, V = 30, CL = 1.7)
+  regim <- new_regimen(amt = sujdos, times = c(0, 12), type= "bolus")
+  out <- sim_ode(ode = mod_1cmt_oral, parameters=param, regimen=regim, t_obs = xtim, only_obs = TRUE)
+  expect_equal(out$t, xtim)
+})
+
+test_that("model size is appropriate (bug: JeHi 20151204)", {
+  pk3cmt <- new_ode_model(
+    code = "
+      dAdt[1] = -KA*A[1];
+      dAdt[2] = KA*A[1] -(Q/V)*A[2] + (Q/V2)*A[3] -(CL/V)*A[2];
+      dAdt[3] = -(Q/V2)*A[3] + (Q/V)*A[2];
+      ",
+    obs = list(cmt = 2, scale = "V")
+  )
+  expect_equal( attr(pk3cmt, "size"), 3)
+})
+
+# ---- Dose Compartment Tests ----
+test_that("Dose is added to correct compartment: specified by model", {
+  set.seed(90)
+  p <- list(CL = 1, V  = 10, KA = 0.5, S2 = .1)
+  r <- new_regimen(amt = 100, times = c(0), type = "infusion")
+  dat <- sim_ode(
+    ode = mod_dose_cmt_2,
+    n_ind = 1,
+    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
+    parameters = p,
+    regimen = r,
+    verbose = FALSE,
+    t_max = 48
+  )
+  # Dose should be in cmt 2
+  expect_equal(dat$y[dat$comp == 1], rep(0, 50))
+  expect_true(all(dat$y[dat$comp == 2][-1] > 0))
+})
+
+test_that("Dose is added to correct compartment: override model by regimen", {
+  set.seed(60)
+  p <- list(CL = 1, V  = 10, KA = 0.5, S2 = .1)
+  r <- new_regimen(
+    amt = c(100, 100, 100),
+    times = c(0, 6, 12),
+    cmt = c(1,2,3),
+    type = "bolus"
+  )
+  dat <- sim_ode(
+    ode = mod_dose_cmt_2,
+    n_ind = 1,
+    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
+    parameters = p,
+    regimen = r,
+    verbose = FALSE,
+    t_max = 48
+  )
+  # Dose should be in cmt 1, 2 and 3
+  expect_true(all(dat$y[dat$comp == 1 & dat$t > 0] > 0))
+  expect_true(max(diff(dat$y[dat$comp == 2])) > 95)
+  expect_true(max(diff(dat$y[dat$comp == 3])) > 95)
+})
+
+test_that("Infusion works for all compartments", {
+  set.seed(44)
+  # Part 1: Specify cmt only with model
+  p <- list(CL = 1, V  = 10, KA = 0.5, S2 = .1)
+  r <- new_regimen(
+    amt = c(100, 100, 100),
+    times = c(0, 6, 12),
+    cmt = c(1,2,3),
+    t_inf = 3,
+    type = "infusion"
+  )
+  dat <- sim_ode(
+    ode = mod_dose_cmt_2,
+    n_ind = 1,
+    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
+    parameters = p,
+    regimen = r,
+    verbose = FALSE,
+    t_max = 48
+  )
+  expect_true(all(dat$y[dat$comp == 1 & dat$t > 0 ] > 0))
+  expect_true(max(diff(dat$y[dat$comp == 2])) > 25)
+  expect_true(max(diff(dat$y[dat$comp == 3])) > 25)
+  expect_equal(round(max(dat$y[dat$comp == 2]), 1), 131.2)
+  expect_equal(round(max(dat$y[dat$comp == 3]), 1), 148.4)
+})
+
+test_that("Duplicate obs returned when specified in arg", {
+  # for first 2 doses, infusion time will just be ignored, but a value has to be specified in the vector
+  p <- list(CL = 1, V  = 10, KA = 0.5, S2=.1)
+  r <- new_regimen(
+    amt = c(100, 100, 100, 100),
+    times = c(0, 6, 12, 18),
+    cmt = c(2, 2, 1, 1),
+    t_inf = c(1, 1, 1, 1),
+    type = c("bolus", "bolus", "infusion", "infusion")
+  )
+  dat <- sim_ode(
+    ode = mod_1cmt_oral,
+    n_ind = 1,
+    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
+    parameters = p,
+    regimen = r,
+    t_obs = c(1, 2, 3, 4, 4, 4, 6), ## see duplicate obs here
+    duplicate_t_obs = T,
+    only_obs = FALSE
+  )
+  expect_equal(length(dat[dat$t == 4,]$y),  9)
+  expect_equal(length(dat$y), 21)
+  expect_equal(sum(is.na(dat$y)), 0)
+})
+
+test_that("Custom t_obs is returned", {
+  t_obs <- seq(from = 0, to = 24, by = .1)
+  p <- list(CL = 1, V  = 10, KA = 0.5, S2=.1)
+  r <- new_regimen(
+    amt = c(100, 100, 100, 100),
+    times = c(0, 6, 12, 18),
+    cmt = c(2, 2, 1, 1),
+    t_inf = c(1, 1, 1, 1),
+    type = c("bolus", "bolus", "infusion", "infusion")
+  )
+  dat <- sim_ode(
+    ode = mod_1cmt_oral,
+    n_ind = 1,
+    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
+    parameters = p,
+    regimen = r,
+    t_obs = t_obs,
+    only_obs = T
+  )
+  expect_equal(mean(diff(t_obs)), mean(diff(dat$t)))
+})
+
+# ---- Covariate Timing Tests ----
+test_that("if covariate time is at end of infusion, end of infusion is still recorded", {
+  # Bug reported by JF
+  pop_est <- list(CL = 1.08, V = 0.98)
+  regimen <- new_regimen(
+    amt = c(1500, 1000, 1500, 1500, 1500, 1500, 1500),
+    type = "infusion",
+    t_inf = c(2, 1, 2, 2, 1, 1, 1),
+    times = c(0, 10.8666666666667, 20.4333333333333, 32.0666666666667, 46.9, 54.9, 62.9 )
+  )
+  covs <- list(
+    WT = new_covariate(value = c(60, 65), times = c(0, 47.9)),
+    CRCL = new_covariate(8), CVVH = new_covariate(0)
+  )
+  pksim <- sim(
+    ode = mod_1cmt_iv,
+    parameters = pop_est,
+    covariates = covs,
+    regimen = regimen,
+    checks = TRUE,
+    only_obs = TRUE
+  )
+  expect_true(all(pksim$y < 1000))
+})
+
+test_that("Covariate table simulation runs", {
+  # this test used to be in the covariate_table_to_list file but
+  # makes more sense here.
+  p <- list(CL = 5, V = 50)
+  reg <- new_regimen (amt = 100, n = 4, interval = 12, type = "bolus",  cmt=1)
+  om <- c(0.01, 1, 0.01)
+  cov_table <- data.frame(
+    id=c(1, 1, 2, 3),
+    WT = c(40, 45, 50, 60),
+    SCR = c(50, 150, 90,110),
+    t = c(0, 5, 0, 0)
+  )
+
+  dat <- sim(
+    mod_1cmt_iv,
+    parameters = p,
+    regimen = reg,
+    covariates_table = cov_table,
+    covariates_implementation = list(SCR = "interpolate"),
+    omega = NULL,
+    n_ind = 3,
+    only_obs = T,
+    output_include = list(parameters = TRUE, covariates=TRUE)
+  )
+  expect_equal(length(unique(dat$id)), 3)
+})
diff --git a/tests/testthat/test_advan.R b/tests/testthat/test_advan.R
deleted file mode 100644
index 5a8dce3a..00000000
--- a/tests/testthat/test_advan.R
+++ /dev/null
@@ -1,126 +0,0 @@
-## These models are also tested in the unit tests for `calc_ss_analytics()`, so just testing a few example cases here
-dose <- 100
-interval <- 12
-t_inf <- 1
-n_days <- 5
-parameters <- list(CL = 10, V = 50, KA = 0.5, Q = 5, V2 = 100, Q2 = 3, V3 = 150, F1 = 1)
-t_obs <- c(3, 6, 8, 23)
-reg_bolus <- new_regimen(
-  amt = dose,
-  times = seq(0, interval * n_days * (24/interval), interval),
-  t_inf = t_inf, type = "bolus"
-)
-data <- advan_create_data(
-  reg_bolus,
-  parameters = parameters,
-  cmts = 5,
-  t_obs = t_obs
-)
-
-## Infusion dataset
-reg_infusion <- new_regimen(
-  amt = dose,
-  times = seq(0, interval * n_days * (24/interval), interval),
-  t_inf = t_inf,
-  type = "infusion"
-)
-data_infusion <- advan_create_data(
-  reg_infusion,
-  parameters = parameters,
-  cmts = 6,
-  t_obs = t_obs
-)
-
-## One compartment
-test_that("One compartment IV bolus", {
-  res1_iv   <- advan("1cmt_iv_bolus", cpp=FALSE)(data)
-  res1_iv_c <- advan("1cmt_iv_bolus", cpp=TRUE)(data)
-  expect_equal(round(res1_iv[res1_iv$TIME == 23,]$DV, 3), 0.242)
-  expect_true(!any(is.na(res1_iv$DV)))
-  expect_equal(res1_iv, res1_iv_c)
-})
-
-test_that("One compartment IV infusion", {
-  res1_iv_inf   <- advan("1cmt_iv_infusion", cpp=FALSE)(data_infusion)
-  res1_iv_inf_c <- advan("1cmt_iv_infusion", cpp=TRUE)(data_infusion)
-  f1 <- advan("1cmt_iv_infusion", cpp=FALSE)
-  f2 <- advan("1cmt_iv_infusion", cpp=TRUE)
-
-  expect_equal(round(res1_iv_inf[res1_iv_inf$TIME == 23,]$DV, 3), 0.268)
-  expect_true(!any(is.na(res1_iv_inf$DV)))
-  expect_equal(res1_iv_inf, res1_iv_inf_c)
-  expect_equal(attr(f1, "type"), "infusion")
-  expect_equal(attr(f2, "type"), "infusion")
-  expect_equal(attr(f1, "implementation"), FALSE)
-  expect_equal(attr(f2, "implementation"), TRUE)
-  expect_equal(attr(f1, "cmt"), 1)
-  expect_equal(attr(f2, "cmt"), 1)
-})
-
-test_that("One compartment oral", {
-  res1_oral <- advan("1cmt_oral", cpp=FALSE)(data)
-  res1_oral_c <- advan("1cmt_oral", cpp=TRUE)(data)
-
-  expect_equal(round(res1_oral[res1_oral$TIME == 23,]$DV, 3), 0.389)
-  expect_true(!any(is.na(res1_oral$DV)))
-  expect_equal(res1_oral, res1_oral_c)
-})
-
-
-## Two compartment
-test_that("Two compartment iv bolus", {
-  res2_iv   <- advan("2cmt_iv_bolus", cpp=FALSE)(data)
-  res2_iv_c <- advan("2cmt_iv_bolus", cpp=TRUE)(data)
-
-  expect_equal(round(res2_iv[res2_iv$TIME == 23,]$DV, 3), 0.212)
-  expect_true(!any(is.na(res2_iv$DV)))
-  expect_equal(res2_iv, res2_iv_c)
-})
-
-test_that("Two compartment iv infusion", {
-  res2_iv_inf   <- advan("2cmt_iv_infusion", cpp=FALSE)(data_infusion)
-  res2_iv_inf_c <- advan("2cmt_iv_infusion", cpp=TRUE)(data_infusion)
-
-  expect_equal(round(res2_iv_inf[res2_iv_inf$TIME == 23,]$DV, 3), 0.225)
-  expect_true(!any(is.na(res2_iv_inf$DV)))
-  expect_equal(res2_iv_inf, res2_iv_inf_c)
-})
-
-test_that("Two compartment oral", {
-  res2_oral   <- advan("2cmt_oral", cpp=FALSE)(data)
-  res2_oral_c <- advan("2cmt_oral", cpp=TRUE)(data)
-
-  expect_equal(round(res2_oral[res2_oral$TIME == 23,]$DV, 3), 0.302)
-  expect_true(!any(is.na(res2_oral$DV)))
-  expect_equal(res2_oral, res2_oral_c)
-})
-
-
-
-## Three compartment
-test_that("Three compartment IV bolus", {
-  res3_iv   <- advan("3cmt_iv_bolus", cpp=FALSE)(data)
-  res3_iv_c <- advan("3cmt_iv_bolus", cpp=TRUE)(data)
-
-  expect_equal(round(res3_iv[res3_iv$TIME == 23,]$DV, 3), 0.169)
-  expect_true(!any(is.na(res3_iv$DV)))
-  expect_equal(res3_iv, res3_iv_c)
-})
-
-test_that("Three compartment IV infusion", {
-  res3_iv_inf   <- advan("3cmt_iv_infusion", cpp=FALSE)(data_infusion)
-  res3_iv_inf_c <- advan("3cmt_iv_infusion", cpp=TRUE)(data_infusion)
-
-  expect_equal(round(res3_iv_inf[res3_iv_inf$TIME == 23,]$DV, 3), 0.177)
-  expect_true(!any(is.na(res3_iv_inf$DV)))
-  expect_equal(res3_iv_inf, res3_iv_inf_c)
-})
-
-test_that("Three compartment IV oral", {
-  res3_oral <- advan("3cmt_oral", cpp=FALSE)(data)
-  res3_oral_c <- advan("3cmt_oral", cpp=TRUE)(data)
-
-  expect_equal(round(res3_oral[res3_oral$TIME == 23,]$DV, 3), 0.236)
-  expect_true(!any(is.na(res3_oral$DV)))
-  expect_equal(res3_oral, res3_oral_c)
-})
diff --git a/tests/testthat/test_advan_with_auc.R b/tests/testthat/test_advan_with_auc.R
deleted file mode 100644
index ced7ef98..00000000
--- a/tests/testthat/test_advan_with_auc.R
+++ /dev/null
@@ -1,173 +0,0 @@
-if (identical(Sys.getenv("NOT_CRAN"), "true")) {
-  dose <- 100
-  interval <- 12
-  n_days <- 5
-  t_inf <- 1.5
-  parameters <- list(
-    CL = 10,
-    V = 50,
-    KA = 0.5,
-    Q = 5,
-    V2 = 100,
-    Q2 = 3,
-    V3 = 150,
-    F1 = 1
-  )
-  t_obs <- c(3, 6, 8, 23, 47)
-
-  ## ODE models for testing
-  mod_1cmt <- new_ode_model(
-    code="dAdt[1] = -(CL/V)*A[1]; dAdt[2] = A[1]/V;",
-    parameters = parameters
-  )
-  mod_2cmt <- new_ode_model(
-    code="
-    dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2];
-    dAdt[2] = +(Q/V)*A[1] - (Q/V2)*A[2];
-    dAdt[3] = A[1]/V;
-  ",
-  parameters = parameters
-  )
-  mod_3cmt <- new_ode_model(
-    code="
-    dAdt[1] = -(CL/V)*A[1] - (Q/V)*A[1] + (Q/V2)*A[2] - (Q2/V)*A[1] + (Q2/V3)*A[3];
-    dAdt[2] =                (Q/V)*A[1]  -(Q/V2)*A[2]                             ;
-    dAdt[3] =                                           (Q2/V)*A[1] - (Q2/V3)*A[3];
-    dAdt[4] = A[1]/V;
-  ",
-  parameters = parameters
-  )
-
-  ## bolus dataset
-  reg_bolus <- new_regimen(
-    amt = dose,
-    times = seq(0, interval * n_days * (24/interval), interval),
-    t_inf = t_inf,
-    type = "bolus"
-  )
-  data_bolus <- advan_create_data(
-    reg_bolus,
-    parameters = parameters,
-    cmts = 5,
-    t_obs = t_obs
-  )
-
-  ## Infusion dataset
-  reg_infusion <- new_regimen(
-    amt = dose,
-    times = seq(0, interval * n_days * (24/interval), interval),
-    t_inf = t_inf,
-    type = "infusion"
-  )
-  data_infusion <- advan_create_data(
-    reg_infusion,
-    parameters = parameters,
-    cmts = 6,
-    t_obs = t_obs
-  )
-}
-
-test_that("One compartment bolus ADVAN runs", {
-  skip_on_cran()
-  res1_iv_r <- advan("1cmt_iv_bolus", cpp=FALSE)(data_bolus)
-  res1_iv_c <- advan("1cmt_iv_bolus", cpp=TRUE)(data_bolus)
-  res1_iv_ode <- sim(ode = mod_1cmt, regimen = reg_bolus, parameters = parameters, t_obs = t_obs)
-
-  # AUC R
-  expect_equal(round(res1_iv_r[res1_iv_r$TIME %in% t_obs,]$AUC, 5), round(res1_iv_ode[res1_iv_ode$comp == 2,]$y, 5))
-
-  #AUC-C
-  expect_equal(round(res1_iv_c[res1_iv_c$TIME %in% t_obs,]$AUC, 5), round(res1_iv_ode[res1_iv_ode$comp == 2,]$y, 5))
-})
-
-test_that("Two compartment bolus ADVAN runs", {
-  skip_on_cran()
-  res2_iv_r <- advan("2cmt_iv_bolus", cpp=FALSE)(data_bolus)
-  res2_iv_c <- advan("2cmt_iv_bolus", cpp=TRUE)(data_bolus)
-  res2_iv_ode <- sim(ode = mod_2cmt, regimen = reg_bolus, parameters = parameters, t_obs = t_obs)
-  expect_equal(
-    round(res2_iv_r[res2_iv_r$TIME %in% t_obs,]$AUC, 5),
-    round(res2_iv_c[res2_iv_c$TIME %in% t_obs,]$AUC, 5)
-  )
-  # AUC R
-  expect_equal(
-    round(res2_iv_r[res2_iv_r$TIME %in% t_obs,]$AUC, 5),
-    round(res2_iv_ode[res2_iv_ode$comp == 3,]$y, 5)
-  )
-
-  #AUC-C
-  expect_equal(
-    round(res2_iv_c[res2_iv_c$TIME %in% t_obs,]$AUC, 5),
-    round(res2_iv_ode[res2_iv_ode$comp == 3,]$y, 5)
-  )
-})
-
-test_that("Two compartment infusion ADVAN runs", {
-  skip_on_cran()
-  res2_inf_r <- advan("2cmt_iv_infusion", cpp=FALSE)(data_infusion)
-  res2_inf_c <- advan("2cmt_iv_infusion", cpp=TRUE)(data_infusion)
-  res2_inf_ode <- sim(ode = mod_2cmt, regimen = reg_infusion, parameters = parameters, t_obs = t_obs)
-
-  expect_equal(
-    round(res2_inf_r[res2_inf_r$TIME %in% t_obs,]$AUC, 5),
-    round(res2_inf_c[res2_inf_c$TIME %in% t_obs,]$AUC, 5)
-  )
-
-  # AUC R
-  expect_equal(
-    round(res2_inf_r[res2_inf_r$TIME %in% t_obs,]$AUC, 5),
-    round(res2_inf_ode[res2_inf_ode$comp == 3,]$y, 5)
-  )
-
-  #AUC-C
-  expect_equal(
-    round(res2_inf_c[res2_inf_c$TIME %in% t_obs,]$AUC, 5),
-    round(res2_inf_ode[res2_inf_ode$comp == 3,]$y, 5)
-  )
-
-})
-
-test_that("Three compartment bolus ADVAN runs", {
-  skip_on_cran()
-  res3_iv_r <- advan("3cmt_iv_bolus", cpp=FALSE)(data_bolus)
-  res3_iv_c <- advan("3cmt_iv_bolus", cpp=TRUE)(data_bolus)
-  res3_iv_ode <- sim(ode = mod_3cmt, regimen = reg_bolus, parameters = parameters, t_obs = t_obs)
-  expect_equal(
-    round(res3_iv_r[res3_iv_r$TIME %in% t_obs,]$AUC, 5),
-    round(res3_iv_c[res3_iv_c$TIME %in% t_obs,]$AUC, 5)
-  )
-  # AUC R
-  expect_equal(
-    round(res3_iv_r[res3_iv_r$TIME %in% t_obs,]$AUC, 5),
-    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
-  )
-
-  #AUC-C
-  expect_equal(
-    round(res3_iv_c[res3_iv_c$TIME %in% t_obs,]$AUC, 5),
-    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
-  )
-})
-
-test_that("Three compartment iv ADVAN runs", {
-  skip_on_cran()
-  res3_iv_r <- advan("3cmt_iv_infusion", cpp=FALSE)(data_infusion)
-  res3_iv_c <- advan("3cmt_iv_infusion", cpp=TRUE)(data_infusion)
-  res3_iv_ode <- sim(ode = mod_3cmt, regimen = reg_infusion, parameters = parameters, t_obs = t_obs)
-  expect_equal(
-    round(res3_iv_r[res3_iv_r$TIME %in% t_obs,]$AUC, 5),
-    round(res3_iv_c[res3_iv_c$TIME %in% t_obs,]$AUC, 5)
-  )
-  # AUC R
-  expect_equal(
-    round(res3_iv_r[res3_iv_r$TIME %in% t_obs,]$AUC, 5),
-    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
-  )
-
-  #AUC-C
-  expect_equal(
-    round(res3_iv_c[res3_iv_c$TIME %in% t_obs,]$AUC, 5),
-    round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
-  )
-})
-
diff --git a/tests/testthat/test_advan_with_covariates.R b/tests/testthat/test_advan_with_covariates.R
deleted file mode 100644
index 7f838358..00000000
--- a/tests/testthat/test_advan_with_covariates.R
+++ /dev/null
@@ -1,62 +0,0 @@
-test_that("Analytic and ODE models with covariates are the same", {
-  skip_on_cran()
-
-  ## Create dataset
-  dose <- 100
-  interval <- 12
-  n_days <- 2
-  parameters <- list(
-    CL = 10,
-    V = 50,
-    KA = 0.5,
-    Q = 5,
-    V2 = 100,
-    Q2 = 3,
-    V3 = 150,
-    F1 = 1
-  )
-  t_obs <- seq(0, 40, .1)
-  reg_bolus <- new_regimen(
-    amt = dose,
-    times = seq(0, interval * n_days * (24/interval), interval),
-    type = "bolus"
-  )
-  ## there is slight difference in how bolus doses are handled.
-  ## Analytical equation is perhaps more consistent, so not testing
-  ## simulations at dose times. Should look into later.
-  t_obs <- t_obs[! t_obs %in% reg_bolus$dose_times]
-  covariates <- list(WT = new_covariate(80), CRCL=new_covariate(4.5))
-
-  ## Using analytic equations model:
-  data_ana <- sim(
-    analytical = "1cmt_iv_bolus",
-    parameters = parameters,
-    covariates = covariates,
-    regimen = reg_bolus,
-    t_obs = t_obs,
-    covariate_model = "CL = CL * (CRCL / 3)^0.75; V = V * (WT / 70.0)"
-  )
-
-  ## Using ODE model:
-  mod1 <- new_ode_model(
-    code = "
-    dAdt[1] = -( (CL*pow(CRCL/3.0, 0.75)) / (V*WT/70.0) ) * A[1];
-  ",
-  covariates = covariates,
-  obs = list(cmt = 1, scale = "V*WT/70.0"), dose = list(cmt = 1)
-  )
-  data_ode <- sim(
-    ode = mod1,
-    parameters = parameters,
-    covariates = covariates,
-    regimen = reg_bolus,
-    t_obs = t_obs,
-    duplicate_t_obs = TRUE,
-    only_obs = TRUE
-  )
-
-  expect_equal(nrow(data_ana), nrow(data_ode))
-  expect_equal(round(data_ana$y,4), round(data_ode$y, 4))
-})
-
-
diff --git a/tests/testthat/test_compare_results.R b/tests/testthat/test_compare_results.R
deleted file mode 100644
index 73652e9f..00000000
--- a/tests/testthat/test_compare_results.R
+++ /dev/null
@@ -1,266 +0,0 @@
-## models: shared between tests and take a while to compile
-#  - oral models
-## Uses model defined in setup.R
-pk1cmt_oral_anal = function(t, dose, KA, V, CL) {
-  dose*KA/(V*(KA-CL/V))*(exp(-(CL/V) * t)-exp(-KA * t))
-}
-pk1cmt_oral_code <- new_ode_model(
-  code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
-  obs=list(cmt = 2, scale="V")
-)
-
-#   - iv models
-## Uses model defined in setup.R
-
-#   - model with dose cmt specified
-dose_in_cmt_2  <- new_ode_model(
-  code = "
-      dAdt[1] = -KA * A[1];
-      dAdt[2] = KA*A[1] -(CL/V) * A[2]
-      dAdt[3] = S2*(A[2]-A[3])
-    ",
-  obs = list(cmt=2, scale="V"),
-  dose = list(cmt = 2),
-  cpp_show_code = FALSE
-)
-
-
-test_that("Library and custom C++ and code matches analytic soln", {
-  p <- list(KA = 1, CL = 5, V = 50)
-  t_obs <- c(0:72)
-  t_obs2 <- t_obs + 0.1234 # also needs to be producing results with non-integer times
-  dose <- 100
-  t_dose <- c(0)
-  regimen <- new_regimen(amt=dose, times = t_dose, type = "oral")
-
-  pk1cmt_oral_lib <- sim_ode(
-    ode = mod_1cmt_oral,
-    parameters = p,
-    regimen = regimen,
-    t_obs = t_obs,
-    int_step_size = 0.1,
-    duplicate_t_obs = TRUE,
-    only_obs=TRUE
-  )
-
-  pk1cmt_oral_code <- sim_ode(
-    ode = pk1cmt_oral_code,
-    parameters = p,
-    duplicate_t_obs = TRUE,
-    regimen=regimen,
-    t_obs=t_obs,
-    int_step_size = 0.1,
-    only_obs=TRUE
-  )
-
-  pk1cmt_oral_anal_res <- pk1cmt_oral_anal(t_obs, dose, p$KA, p$V, p$CL)
-  expect_equal(round(pk1cmt_oral_lib$y, 3), round(pk1cmt_oral_anal_res, 3))
-  expect_equal(round(pk1cmt_oral_code$y, 3), round(pk1cmt_oral_anal_res, 3))
-})
-
-
-test_that("precision in time does not impact # obs returned", {
-  regimen_mult <- new_regimen(
-    amt = rep(12.8, 3),
-    times = c(0, 6, 12),
-    type="infusion",
-    t_inf = 2
-  )
-  t_obs <- c(11.916, 14.000, 16.000, 17.000, 30)
-  tmp <- sim_ode(
-    ode = mod_1cmt_iv,
-    parameters = list(CL = 5, V = 50),
-    regimen = regimen_mult,
-    t_obs = t_obs,
-    only_obs = TRUE
-  )
-  expect_equal(tmp$t, t_obs)
-})
-
-test_that("test bug EmCo 20150925", {
-  xtim <- c(0, 2, 4, 8, 12, 24)
-  sujdos <- 320
-  param <- list(KA = 1.8, V = 30, CL = 1.7)
-  regim <- new_regimen(amt = sujdos, times = c(0, 12), type= "bolus")
-  out <- sim_ode(ode = mod_1cmt_oral, parameters=param, regimen=regim, t_obs = xtim, only_obs = TRUE)
-  expect_equal(out$t, xtim)
-})
-
-test_that("model size is appropriate (bug: JeHi 20151204)", {
-  pk3cmt <- new_ode_model(
-    code = "
-      dAdt[1] = -KA*A[1];
-      dAdt[2] = KA*A[1] -(Q/V)*A[2] + (Q/V2)*A[3] -(CL/V)*A[2];
-      dAdt[3] = -(Q/V2)*A[3] + (Q/V)*A[2];
-      ",
-    obs = list(cmt = 2, scale = "V")
-  )
-  expect_equal( attr(pk3cmt, "size"), 3)
-})
-
-test_that("Dose is added to correct compartment: specified by model", {
-  set.seed(90)
-  p <- list(CL = 1, V  = 10, KA = 0.5, S2 = .1)
-  r <- new_regimen(amt = 100, times = c(0), type = "infusion")
-  dat <- sim_ode(
-    ode = dose_in_cmt_2,
-    n_ind = 1,
-    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
-    parameters = p,
-    regimen = r,
-    verbose = FALSE,
-    t_max = 48
-  )
-  # Dose should be in cmt 2
-  expect_equal(dat$y[dat$comp == 1], rep(0, 50))
-  expect_true(all(dat$y[dat$comp == 2][-1] > 0))
-})
-
-test_that("Dose is added to correct compartment: override model by regimen", {
-  set.seed(60)
-  p <- list(CL = 1, V  = 10, KA = 0.5, S2 = .1)
-  r <- new_regimen(
-    amt = c(100, 100, 100),
-    times = c(0, 6, 12),
-    cmt = c(1,2,3),
-    type = "bolus"
-  )
-  dat <- sim_ode(
-    ode = dose_in_cmt_2,
-    n_ind = 1,
-    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
-    parameters = p,
-    regimen = r,
-    verbose = FALSE,
-    t_max = 48
-  )
-  # Dose should be in cmt 1, 2 and 3
-  expect_true(all(dat$y[dat$comp == 1 & dat$t > 0] > 0))
-  expect_true(max(diff(dat$y[dat$comp == 2])) > 95)
-  expect_true(max(diff(dat$y[dat$comp == 3])) > 95)
-})
-
-test_that("Infusion works for all compartments", {
-  set.seed(44)
-  # Part 1: Specify cmt only with model
-  p <- list(CL = 1, V  = 10, KA = 0.5, S2 = .1)
-  r <- new_regimen(
-    amt = c(100, 100, 100),
-    times = c(0, 6, 12),
-    cmt = c(1,2,3),
-    t_inf = 3,
-    type = "infusion"
-  )
-  dat <- sim_ode(
-    ode = dose_in_cmt_2,
-    n_ind = 1,
-    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
-    parameters = p,
-    regimen = r,
-    verbose = FALSE,
-    t_max = 48
-  )
-  expect_true(all(dat$y[dat$comp == 1 & dat$t > 0 ] > 0))
-  expect_true(max(diff(dat$y[dat$comp == 2])) > 25)
-  expect_true(max(diff(dat$y[dat$comp == 3])) > 25)
-  expect_equal(round(max(dat$y[dat$comp == 2]), 1), 131.2)
-  expect_equal(round(max(dat$y[dat$comp == 3]), 1), 148.4)
-})
-
-test_that("Duplicate obs returned when specified in arg", {
-  # for first 2 doses, infusion time will just be ignored, but a value has to be specified in the vector
-  p <- list(CL = 1, V  = 10, KA = 0.5, S2=.1)
-  r <- new_regimen(
-    amt = c(100, 100, 100, 100),
-    times = c(0, 6, 12, 18),
-    cmt = c(2, 2, 1, 1),
-    t_inf = c(1, 1, 1, 1),
-    type = c("bolus", "bolus", "infusion", "infusion")
-  )
-  dat <- sim_ode(
-    ode = mod_1cmt_oral,
-    n_ind = 1,
-    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
-    parameters = p,
-    regimen = r,
-    t_obs = c(1, 2, 3, 4, 4, 4, 6), ## see duplicate obs here
-    duplicate_t_obs = T,
-    only_obs = FALSE
-  )
-  expect_equal(length(dat[dat$t == 4,]$y),  9)
-  expect_equal(length(dat$y), 21)
-  expect_equal(sum(is.na(dat$y)), 0)
-})
-
-test_that("Custom t_obs is returned", {
-  t_obs <- seq(from = 0, to = 24, by = .1)
-  p <- list(CL = 1, V  = 10, KA = 0.5, S2=.1)
-  r <- new_regimen(
-    amt = c(100, 100, 100, 100),
-    times = c(0, 6, 12, 18),
-    cmt = c(2, 2, 1, 1),
-    t_inf = c(1, 1, 1, 1),
-    type = c("bolus", "bolus", "infusion", "infusion")
-  )
-  dat <- sim_ode(
-    ode = mod_1cmt_oral,
-    n_ind = 1,
-    omega = cv_to_omega(par_cv = list("CL"=0.1, "V"=0.1, "KA" = .1), p),
-    parameters = p,
-    regimen = r,
-    t_obs = t_obs,
-    only_obs = T
-  )
-  expect_equal(mean(diff(t_obs)), mean(diff(dat$t)))
-})
-
-test_that("if covariate time is at end of infusion, end of infusion is still recorded", {
-  # Bug reported by JF
-  pop_est <- list(CL = 1.08, V = 0.98)
-  regimen <- new_regimen(
-    amt = c(1500, 1000, 1500, 1500, 1500, 1500, 1500),
-    type = "infusion",
-    t_inf = c(2, 1, 2, 2, 1, 1, 1),
-    times = c(0, 10.8666666666667, 20.4333333333333, 32.0666666666667, 46.9, 54.9, 62.9 )
-  )
-  covs <- list(
-    WT = new_covariate(value = c(60, 65), times = c(0, 47.9)),
-    CRCL = new_covariate(8), CVVH = new_covariate(0)
-  )
-  pksim <- sim(
-    ode = mod_1cmt_iv,
-    parameters = pop_est,
-    covariates = covs,
-    regimen = regimen,
-    checks = TRUE,
-    only_obs = TRUE
-  )
-  expect_true(all(pksim$y < 1000))
-})
-
-test_that("Covariate table simulation runs", {
-  # this test used to be in the covariate_table_to_list file but
-  # makes more sense here.
-  p <- list(CL = 5, V = 50)
-  reg <- new_regimen (amt = 100, n = 4, interval = 12, type = "bolus",  cmt=1)
-  om <- c(0.01, 1, 0.01)
-  cov_table <- data.frame(
-    id=c(1, 1, 2, 3),
-    WT = c(40, 45, 50, 60),
-    SCR = c(50, 150, 90,110),
-    t = c(0, 5, 0, 0)
-  )
-
-  dat <- sim(
-    mod_1cmt_iv,
-    parameters = p,
-    regimen = reg,
-    covariates_table = cov_table,
-    covariates_implementation = list(SCR = "interpolate"),
-    omega = NULL,
-    n_ind = 3,
-    only_obs = T,
-    output_include = list(parameters = TRUE, covariates=TRUE)
-  )
-  expect_equal(length(unique(dat$id)), 3)
-})
diff --git a/tests/testthat/test_iov.R b/tests/testthat/test_iov.R
index 650d1388..eb757c00 100644
--- a/tests/testthat/test_iov.R
+++ b/tests/testthat/test_iov.R
@@ -1,3 +1,4 @@
+# Uses mod_1cmt_no_iov and mod_1cmt_iov from setup.R
 pars <- list(
   "kappa_CL_1" = 0,
   "kappa_CL_2" = 0,
@@ -12,32 +13,6 @@ pars0 <- list(
   "V" = 50,
   "KA" = 1
 )
-pk0 <- new_ode_model( # no IOV
-  code = "
-      dAdt[1] = -KA * A[1]
-      dAdt[2] = +KA * A[1] -(CL/V) * A[2]
-    ",
-  obs = list(cmt = 2, scale = "V"),
-  dose = list(cmt = 1, bioav = 1),
-  parameters = names(pars0),
-  cpp_show_code = F
-)
-pk1 <- new_ode_model(
-  code = "
-      CL_iov = CL * exp(kappa_CL + eta_CL);
-      dAdt[1] = -KA * A[1]
-      dAdt[2] = +KA * A[1] -(CL_iov/V) * A[2]
-    ",
-  iov = list(
-    cv = list(CL = 0.2),
-    n_bins = 3
-  ),
-  obs = list(cmt = 2, scale = "V"),
-  dose = list(cmt = 1, bioav = 1),
-  declare_variables = c("kappa_CL", "CL_iov"),
-  parameters = names(pars),
-  cpp_show_code = F
-)
 reg1 <- new_regimen(
   amt = 100,
   interval = 24,
@@ -49,7 +24,7 @@ iov_var <- 0.3 ^ 2 # 30% IOV
 test_that("Throws error when `iov_bins` supplied but not present in model", {
   expect_error({
     sim(
-      ode = pk0,
+      ode = mod_1cmt_no_iov,
       parameters = pars0,
       regimen = reg1,
       omega = c(
@@ -67,7 +42,7 @@ test_that("Throws error when `iov_bins` supplied but not present in model", {
 test_that("Throws error when number of `iov_bins` is higher than allowed for model", {
   expect_error({
     sim(
-      ode = pk1,
+      ode = mod_1cmt_iov,
       parameters = pars,
       regimen = reg1,
       omega = c(
@@ -89,7 +64,7 @@ test_that("Throws error when number of `iov_bins` is higher than allowed for mod
 test_that("Throws warning when number of `iov_bins` is lower than allowed for model", {
   expect_warning({
     sim(
-      ode = pk1,
+      ode = mod_1cmt_iov,
       parameters = pars,
       regimen = reg1,
       omega = c(
@@ -113,7 +88,7 @@ test_that("IOV is added to parameters", {
   set.seed(32)
 
   dat <- sim(
-    ode = pk1,
+    ode = mod_1cmt_iov,
     parameters = pars,
     regimen = reg1,
     omega = c(
@@ -217,7 +192,7 @@ test_that("error is not invoked when using parameters_table", {
   # specifying both parameters_table but for a model with IOV should not fail!
   expect_silent(
     dat <- sim(
-      ode = pk1,
+      ode = mod_1cmt_iov,
       parameters_table = parameters_table,
       regimen = reg1,
       omega = c(
@@ -238,7 +213,7 @@ test_that("error is not invoked when using parameters_table", {
   # specifying both parameters and parameters_table should fail
   expect_error(
     dat <- sim(
-      ode = pk1,
+      ode = mod_1cmt_iov,
       parameters = pars,
       parameters_table = parameters_table,
       regimen = reg1,
diff --git a/tests/testthat/test_join_regimen.R b/tests/testthat/test_join_regimen.R
deleted file mode 100644
index 193a2c85..00000000
--- a/tests/testthat/test_join_regimen.R
+++ /dev/null
@@ -1,93 +0,0 @@
-reg1 <- new_regimen(
-  amt = c(100, 100, 100),
-  t_inf = c(24, 24, 28),
-  time = c(0, 24, 48),
-  type = "infusion"
-)
-reg2 <- new_regimen(
-  amt = c(200, 200, 200),
-  t_inf = c(24,24,24),
-  time = c(0, 24, 48),
-  type = "infusion"
-)
-
-test_that("Early update of a regimen works", {
-  res1 <- join_regimen(reg1, reg2, t_dose_update = 65)
-  res2 <- join_regimen(reg1, reg2, interval = 15)
-  expect_equal(res1$t_inf, c(24, 24, 17, 24, 24, 24))
-  expect_equal(round(res1$dose_amts[3], 1), round(100 * 17 / 28, 1))
-  expect_equal(res2$t_inf, c(24, 24, 15, 24, 24, 24))
-  expect_equal(round(res2$dose_amts[3], 1), round(100 * 15 / 28, 1))
-})
-
-test_that("Gap in regimen is ok", {
-  res3 <- join_regimen(reg1, reg2, t_dose_update = 100)
-  expect_equal(res3$t_inf, c(24, 24, 28, 24, 24, 24))
-  expect_equal(res3$dose_amts, c(100, 100, 100, 200, 200, 200))
-})
-
-test_that("Update from time = 0", {
-  res4 <- join_regimen(reg1, reg2, t_dose_update = 0)
-  expect_equal(res4$dose_amts,reg2$dose_amts)
-})
-
-
-test_that("Early update of a regimen works when mixed routes", {
-  reg3 <- new_regimen(
-    amt = c(100, 100, 100),
-    t_inf = c(0, 0, 0),
-    time = c(0, 24, 48),
-    type = "oral"
-  )
-  res0 <- join_regimen(reg3, reg2, dose_update = 1)
-  expect_equal(res0$type, rep("infusion", 3))
-
-  res1 <- join_regimen(reg3, reg2, dose_update = 3)
-  expect_equal(res1$type, c(rep("oral", 2), rep("infusion", 3)))
-
-  res2 <- join_regimen(reg3, reg2, dose_update = 4, interval = 24)
-  expect_equal(res2$type, c(rep("oral", 3), rep("infusion", 3)))
-})
-
-test_that("insufficient join timing throws errors", {
-  expect_error(
-    join_regimen(reg1, reg2)
-  )
-})
-
-test_that("dose_update longer than reg 1 length without interval specified throws error", {
-  expect_error(
-    join_regimen(reg1, reg2, dose_update = 4)
-  )
-})
-
-test_that("one null regimen returns the other", {
-  expect_equal(join_regimen(NULL, reg2), reg2)
-  expect_equal(join_regimen(reg1, NULL), reg1)
-})
-
-test_that("join regimen works when no t_inf specified, e.g. oral or sc dosing", {
-  reg_i <- new_regimen(amt = 500, type = "sc", times = c(0, 2, 6)*24*7)
-  reg_m <- new_regimen(amt = 400, type = "sc", n = 5, interval = 8*7*24)
-  reg <- join_regimen(reg_i, reg_m, interval = 8*7*24)
-  expect_equal(reg$dose_times, c(0, 336, 1008, 2352, 3696, 5040, 6384, 7728))
-  expect_equal(reg$type, c("sc", "sc", "sc", "sc", "sc", "sc", "sc", "sc"))
-  expect_equal(reg$dose_amts, c(500, 500, 500, 400, 400, 400, 400, 400))
-})
-
-test_that("join_regimen with t_dose_update maintains correct rate length", {
-  reg1 <- new_regimen(
-    amt = 1,
-    times = c(0, 24),
-    type = "oral"
-  )
-  reg2 <- new_regimen(
-    amt = 2,
-    time = c(0, 24, 48),
-    type = "oral"
-  )
-  reg3 <- join_regimen(reg1, reg2, t_dose_update = 24)
-
-  expect_equal(length(reg3$rate), length(reg3$dose_times))
-  expect_equal(reg3$n, 4)
-})
diff --git a/tests/testthat/test_merge_regimen.R b/tests/testthat/test_merge_regimen.R
deleted file mode 100644
index fa0bb904..00000000
--- a/tests/testthat/test_merge_regimen.R
+++ /dev/null
@@ -1,32 +0,0 @@
-test_that("merge regimens correctly", {
-  reg1 <- new_regimen(amt = 1000, times = c(0, 24), type = "infusion", t_inf = 1)
-  reg2 <- new_regimen(amt = 500, times = c(12, 36), type = "oral")
-  reg3 <- merge_regimen(regimens = list(reg1, reg2))
-  expect_equal(reg3$dose_times, c(0, 12, 24, 36))
-  expect_equal(reg3$type, c("infusion", "oral", "infusion", "oral"))
-  expect_equal(reg3$dose_amts, c(1000, 500, 1000, 500))
-  expect_equal(reg3$t_inf, c(1, 0, 1, 0))
-  expect_null(reg3$cmt)
-})
-
-test_that("merge regimens correctly: two infusions", {
-  reg1 <- new_regimen(amt = 1000, times = c(0, 24), type = "infusion_1", t_inf = 1)
-  reg2 <- new_regimen(amt = 500, times = c(12, 36), type = "infusion_2", t_inf = 2)
-  reg3 <- merge_regimen(regimens = list(reg1, reg2))
-  expect_equal(reg3$dose_times, c(0, 12, 24, 36))
-  expect_equal(reg3$type, c("infusion_1", "infusion_2", "infusion_1", "infusion_2"))
-  expect_equal(reg3$dose_amts, c(1000, 500, 1000, 500))
-  expect_equal(reg3$t_inf, c(1, 2, 1, 2))
-  expect_null(reg3$cmt)
-})
-
-test_that("cmt info gets merged as well, when available", {
-  reg1 <- new_regimen(amt = 1000, times = c(0, 24), type = "infusion", t_inf = 1, cmt = 2)
-  reg2 <- new_regimen(amt = 500, times = c(12, 36), type = "oral", cmt = 1)
-  reg3 <- merge_regimen(regimens = list(reg1, reg2))
-  expect_equal(reg3$dose_times, c(0, 12, 24, 36))
-  expect_equal(reg3$type, c("infusion", "oral", "infusion", "oral"))
-  expect_equal(reg3$dose_amts, c(1000, 500, 1000, 500))
-  expect_equal(reg3$t_inf, c(1, 0, 1, 0))
-  expect_equal(reg3$cmt, c(2, 1, 2, 1))
-})
diff --git a/tests/testthat/test_new_regimen.R b/tests/testthat/test_new_regimen.R
deleted file mode 100644
index 05a685eb..00000000
--- a/tests/testthat/test_new_regimen.R
+++ /dev/null
@@ -1,123 +0,0 @@
-test_that("Basic regimen creation working", {
-  reg <- new_regimen(amt=100, n=4, interval=4)
-  expect_true(all(c("regimen", "list") %in% class(reg)))
-})
-
-test_that("Regimen type parsed correctly", {
-  reg_oral <- new_regimen(amt = 100, n = 4, interval = 4, type = "oral")
-  reg_bolus <- new_regimen(amt = 100, n=4, interval = 4, type = "bolus")
-  reg_inf <- new_regimen(amt = 100, n = 4, interval = 4, type = "infusion")
-  reg_sc <- new_regimen(amt = 100, n = 4, interval = 4, type = "sc")
-  reg_im <- new_regimen(amt = 100, n = 4, interval = 4, type = "im")
-  reg_mixed <- new_regimen(amt = 100, n = 5, interval = 4, type = c("bolus", "infusion", "oral", "sc", "im"))
-  reg_oral_susp <- new_regimen(amt = 1, n = 4, interval = 8, type = "oral_susp")
-
-  expect_equal(reg_oral$type, rep("oral", 4))
-  expect_equal(reg_bolus$type, rep("bolus", 4))
-  expect_equal(reg_inf$type, rep("infusion", 4))
-  expect_equal(reg_sc$type, rep("sc", 4))
-  expect_equal(reg_im$type, rep("im", 4))
-  expect_equal(reg_mixed$type, c("bolus", "infusion", "oral", "sc", "im"))
-  expect_equal(reg_oral_susp$type, rep("oral_susp", 4))
-
-  # oral-type regimens should have t_inf set to 0
-  expect_equal(reg_oral$t_inf, rep(0, 4))
-  expect_equal(reg_oral_susp$t_inf, rep(0, 4))
-})
-
-test_that("Auto-detect infusion vs bolus", {
-  expect_silent(reg1 <- new_regimen(amt = 100, n = 4, interval = 4, t_inf = 0))
-  expect_warning(reg2 <- new_regimen(amt = 100, n = 4, interval = 4, t_inf = 1))
-  expect_warning(reg3 <- new_regimen(amt = 100, n = 4, interval = 4, t_inf = c(0, 1, 0, 1)))
-
-  expect_equal(reg1$type, rep("bolus", 4))
-  expect_equal(reg2$type, rep("infusion", 4))
-  expect_equal(reg3$type, c("bolus", "infusion", "bolus", "infusion"))
-})
-
-test_that("n = 0 doses does not create doses at negative times", {
-  reg0 <- new_regimen(amt = 200, n = 0, interval = 24, t_inf = 2, type = "infusion")
-  expect_false(min(reg0$dose_times) < 0)
-})
-
-test_that("Rate argument creates valid PKPD regimen", {
-  reg1 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36, 48),
-    type = "infusion",
-    rate = c(1,2,3,4,5)
-  )
-  expect_equal(round(reg1$t_inf), c(100, 50, 33, 25, 20))
-
-  reg2 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36, 48),
-    type = "infusion",
-    rate = c(5)
-  )
-  expect_equal(round(reg2$t_inf), rep(20, 5))
-})
-
-test_that("Doses < 0 set to 0", {
-  expect_warning(
-    tmp <- new_regimen(amt = c(-1, -2, 3, 4), times = c(0, 24, 48, 72), type = "infusion")
-  )
-  expect_true(all(tmp$dose_amts >= 0))
-})
-
-test_that("new_regimen can take arbitrary values for `type`", {
-  reg <- new_regimen(100, times = 0, type = "pip")
-  expect_equal(reg$type, "pip")
-})
-
-test_that("do not creat regimens of `type` 'covariate'", {
-  expect_error(new_regimen(100, times = 0, type = "covariate"))
-})
-
-test_that("sc doses accept an infusion length argument'", {
-  reg1 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36, 48),
-    type = "sc",
-    t_inf = 30/60
-  )
-  expect_equal(reg1$t_inf, rep(0.5,5))
-})
-
-test_that("t_inf imputed correctly", {
-  reg1 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36, 48, 60, 72, 84),
-    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral")
-  )
-  reg2 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36, 48, 60, 72, 84),
-    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral"),
-    t_inf = c(2/60, 2.5, 3/60, NA, NA, NA, NA, NA)
-  )
-  reg3 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36, 48, 60, 72, 84),
-    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral"),
-    t_inf = numeric(0)
-  )
-  reg4 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36, 48, 60, 72, 84),
-    type = c("sc", "infusion", "im", "sc", "infusion", "im","bolus","oral"),
-    t_inf = NULL
-  )
-  reg5 <- new_regimen(
-    amt = 100,
-    times = c(0, 12, 24, 36),
-    type = c("sc", "infusion", "im", "unknown_drug_type"),
-    t_inf = c(2/60, 2.5, 3/60, NA)
-  )
-  expect_equal(reg1$t_inf, c(1/60, 1, 1/60, 1/60, 1, 1/60, 0, 0))
-  expect_equal(reg2$t_inf, c(2/60, 2.5, 3/60, 1/60, 1, 1/60, 0, 0))
-  expect_equal(reg3$t_inf, c(1/60, 1, 1/60, 1/60, 1, 1/60, 0, 0))
-  expect_equal(reg4$t_inf, c(1/60, 1, 1/60, 1/60, 1, 1/60, 0, 0))
-  expect_equal(reg5$t_inf, c(2/60, 2.5, 3/60, 1))
-})
-
diff --git a/tests/testthat/test_shift_regimen.R b/tests/testthat/test_shift_regimen.R
deleted file mode 100644
index 8cee6b56..00000000
--- a/tests/testthat/test_shift_regimen.R
+++ /dev/null
@@ -1,31 +0,0 @@
-test_that("Shifting by 1 works", {
-  keys <- c("oral", "infusion", "bolus")
-  for(key in keys) {
-    reg1 <- new_regimen(amt = 2000, interval = 24, n = 6, type = key, t_inf = 1)
-    reg1_s1 <- shift_regimen(reg1)
-    expect_true("regimen" %in% class(reg1_s1))
-    expect_equal(length(reg1_s1$dose_amts), 5)
-  }
-})
-
-test_that("Shifting by N works", {
-  keys <- c("oral", "infusion", "bolus")
-  for(key in keys) {
-    reg2 <- new_regimen(amt = c(1:6), times = c(0:5) * 24, type = "infusion")
-    reg2_s1 <- shift_regimen(reg2, n = 3)
-    expect_true("regimen" %in% class(reg2_s1))
-    expect_equal(length(reg2_s1$dose_amts), 3)
-    expect_equal(length(reg2_s1$t_inf), 3)
-    # amounts taken from end, not front
-    expect_equal(reg2_s1$dose_amts, c(4, 5, 6))
-  }
-})
-
-test_that("Shifting by N > length(regimen) produces NULL", {
-  keys <- c("oral", "infusion", "bolus")
-  for(key in keys) {
-    reg2 <- new_regimen(amt = c(1:6), times = c(0:5) * 24, type = "infusion")
-    reg2_s2 <- shift_regimen(reg2, n = 10)
-    expect_null(reg2_s2)
-  }
-})
diff --git a/tests/testthat/test_sim_core.R b/tests/testthat/test_sim_core.R
deleted file mode 100644
index fd81ca48..00000000
--- a/tests/testthat/test_sim_core.R
+++ /dev/null
@@ -1,51 +0,0 @@
-test_that("sim core works", {
-  # Uses model defined in setup.R
-
-  reg <- new_regimen(amt = 100, n = 5, interval = 12, t_inf = 1, type = "infusion")
-  par <- list(KA = 1, CL = 5, V = 50)
-
-  ## have to be explicit about t_obs with sim_core!
-  f1 <- function() {
-    res <- sim(ode = mod_1cmt_oral, regimen = reg, parameters = par, only_obs = TRUE, t_obs=c(0:24))$y
-    return(res)
-  }
-
-  f2 <- function() {
-    obj <- sim(ode = mod_1cmt_oral, regimen = reg, parameters = par, only_obs = TRUE, t_obs=c(0:24), return_design=TRUE)
-    sim_core(obj, ode = mod_1cmt_oral)$y
-  }
-  expect_equal(f1(), f2())
-})
-
-test_that("sim core works for absorption model with lagtime", {  
-  reg <- new_regimen(amt = 100, n = 5, interval = 12, t_inf = 1, type = "infusion")
-  par <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
-
-  ## have to be explicit about t_obs with sim_core!
-  f_ref <- function() {
-    res <- sim(
-      ode = mod_1cmt_oral_lagtime, 
-      regimen = reg, 
-      parameters = par, 
-      only_obs = TRUE, 
-      t_obs=c(0:24)
-    )$y
-    return(res)
-  }
-  f_core <- function() {
-    obj <- sim(
-      ode = mod_1cmt_oral_lagtime, 
-      regimen = reg, 
-      parameters = par, 
-      only_obs = TRUE, 
-      t_obs=c(0:24), 
-      return_design=TRUE
-    )
-    sim_core(
-      obj, 
-      ode = mod_1cmt_oral_lagtime, 
-      lagtime = c(0.83, 0)
-    )$y ## Lagtime parameter needed!
-  }
-  expect_equal(f_ref(), f_core())
-})
diff --git a/tests/testthat/test_sim_lagtime.R b/tests/testthat/test_sim_lagtime.R
deleted file mode 100644
index 77a8c5f7..00000000
--- a/tests/testthat/test_sim_lagtime.R
+++ /dev/null
@@ -1,26 +0,0 @@
-test_that("dose dump after lagtime in correct order in output data", {
-  skip_on_cran()
-  reg <- new_regimen(amt = 500, n = 4, interval = 12, type = 'oral')
-  pars <- list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
-  dat <- sim_ode(
-    ode = mod_1cmt_oral_lagtime,
-    regimen = reg,
-    parameters = pars,
-    only_obs = FALSE
-  )
-  ## Change after RXR-2394: time of TLAG not in dataset anymore unless requested by user in t_obs
-  ## Before change: expect_equal(round(dat[dat$t == 12.83 & dat$comp == 1,]$y, 1), c(1.2, 501.2))
-  ## After change:
-  expect_equal(nrow(dat[dat$t == 12.83,]), 0)
-  ## When grid requested by user, lagtime should be visible
-  dat <- sim_ode(
-    ode = mod_1cmt_oral_lagtime,
-    regimen = reg,
-    parameters = pars,
-    t_obs = seq(0, 1, 0.01),
-    only_obs = TRUE
-  )
-  tmp <- dat[dat$t >= 0.82 & dat$t <= 0.85, ]
-  expect_equal(tmp$t, c(0.82, 0.83, 0.84, 0.85))
-  expect_equal(round(tmp$y, 2), c(0, 0, 0.05, 0.10))
-})