From 10fc76d0d37de09538eb225e13cadd32f932fe9d Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Thu, 15 Jan 2026 20:39:43 -0800
Subject: [PATCH 01/10] consolidate model compilations in test setup.R

Move shared model compilations from individual test files to setup.R
to reduce redundant compilation during test runs. Models are compiled
once at test session startup and reused across test files.

Models moved include IOV models, multi-observation models, mixture
models, time-varying covariate models, and various library models
(3cmt, 2cmt oral, MM kinetics). Conditional compilation (NOT_CRAN)
preserved for slower models.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 tests/testthat/setup.R                  | 182 +++++++++++++++++++++++-
 tests/testthat/test_advan_with_auc.R    |  52 ++-----
 tests/testthat/test_calc_ss_analytic.R  |  13 +-
 tests/testthat/test_cmt_mapping.R       |  13 +-
 tests/testthat/test_compare_results.R   |  32 ++---
 tests/testthat/test_get_var_y.R         |  12 +-
 tests/testthat/test_iov.R               | 107 +++++---------
 tests/testthat/test_mixture_model.R     |  37 ++---
 tests/testthat/test_multi_obs.R         |  59 ++++----
 tests/testthat/test_reparametrization.R |  40 ++----
 tests/testthat/test_t_init.R            |  26 ++--
 tests/testthat/test_timevar_cov.R       |  35 ++---
 12 files changed, 330 insertions(+), 278 deletions(-)

diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index 906863f0..ed1138f5 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -1,6 +1,13 @@
+# =============================================================================
+# Shared models compiled once at test session startup
+# =============================================================================
+
+# --- Library models (standard PK 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")
+
+# --- Custom 1cmt oral with lagtime ---
 mod_1cmt_oral_lagtime <- new_ode_model(
   code = "
     dAdt[0] = -KA * A[0]
@@ -11,7 +18,9 @@ mod_1cmt_oral_lagtime <- new_ode_model(
   dose = list(cmt = 1, bioav = 1),
   parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
 )
-oral_1cmt_allometric <- new_ode_model( # also timevarying and dose-dependence factor
+
+# --- 1cmt oral with allometric scaling, time-varying CL, dose-dependent V ---
+oral_1cmt_allometric <- new_ode_model(
   code = "
       if(t<168.0) {
         CLi = CL * pow(WT/70, 0.75)
@@ -32,3 +41,174 @@ oral_1cmt_allometric <- new_ode_model( # also timevarying and dose-dependence fa
   declare_variables = c("CLi", "Vi"),
   parameters = list(KA = 0.5, CL = 5, V = 50)
 )
+
+# --- From test_iov.R: IOV models ---
+pars_iov <- list(
+  "kappa_CL_1" = 0,
+  "kappa_CL_2" = 0,
+  "kappa_CL_3" = 0,
+  "eta_CL" = 0,
+  "CL" = 5,
+  "V" = 50,
+  "KA" = 1
+)
+pars_iov_no_iov <- list(
+  "CL" = 5,
+  "V" = 50,
+  "KA" = 1
+)
+pk_iov_none <- 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 = names(pars_iov_no_iov),
+  cpp_show_code = FALSE
+)
+pk_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 = names(pars_iov),
+  cpp_show_code = FALSE
+)
+
+# --- From test_compare_results.R ---
+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
+)
+
+# --- From test_cmt_mapping.R: 1cmt oral with cmt_mapping ---
+pk1cmt_oral_cmt_mapping <- new_ode_model(
+  code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
+  obs = list(cmt = 2, scale = "V"),
+  cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
+)
+
+# --- From test_multi_obs.R: Multi-observation model ---
+vars_multi_obs <- c("CONC", "METAB", "METAB2", "ACT")
+pk_multi_obs <- new_ode_model(
+  code = "dAdt[1] = -(CL/V)*A[1]; CONC = 1000*A[1]/V; METAB = CONC/2; METAB2 = CONC * t; ACT = 15",
+  obs = list(variable = vars_multi_obs),
+  declare_variables = vars_multi_obs,
+  cpp_show_code = FALSE
+)
+
+# --- From test_mixture_model.R ---
+covs_mixture <- list(WT = new_covariate(70))
+mod_mixture <- new_ode_model(
+  code = "
+      dAdt[0] = -(CL*(WT/70.0)/V)*A[0];
+    ",
+  pk_code = " ",
+  obs = list(cmt = 1, scale = "V"),
+  mixture = list(CL = list(values = c(5, 15), probability = 0.3)),
+  covariates = covs_mixture
+)
+
+# --- Conditional models (skip on CRAN due to compilation time) ---
+if (identical(Sys.getenv("NOT_CRAN"), "true")) {
+
+  # Library models
+  pk_3cmt_iv <- new_ode_model("pk_3cmt_iv")
+  pk_2cmt_oral <- new_ode_model("pk_2cmt_oral")
+  pk_3cmt_oral <- new_ode_model("pk_3cmt_oral")
+  mod_1cmt_iv_mm <- new_ode_model("pk_1cmt_iv_mm")
+
+  # --- From test_advan_with_auc.R: Models with AUC compartments ---
+  parameters_advan_auc <- list(
+    CL = 10, V = 50, KA = 0.5, Q = 5, V2 = 100, Q2 = 3, V3 = 150, F1 = 1
+  )
+  mod_1cmt_auc <- new_ode_model(
+    code = "dAdt[1] = -(CL/V)*A[1]; dAdt[2] = A[1]/V;",
+    parameters = parameters_advan_auc
+  )
+  mod_2cmt_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 = parameters_advan_auc
+  )
+  mod_3cmt_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 = parameters_advan_auc
+  )
+
+  # --- From test_timevar_cov.R: 2cmt with time-varying covariates ---
+  mod_2cmt_timevar <- new_ode_model(
+    code = "
+      dAdt[1] = -(Q/V)*A[1] + (Q/V2)*A[2] -(CLi/V)*A[1];
+      dAdt[2] = -(Q/V2)*A[2] + (Q/V)*A[1];
+    ",
+    pk_code = "CLi = CL + CRCL",
+    obs = list(cmt = 2, scale = "V"),
+    covariates = list(CRCL = new_covariate(5)),
+    declare_variables = "CLi",
+    cpp = FALSE
+  )
+
+  # --- From test_t_init.R: Model with state initialization ---
+  mod_t_init <- new_ode_model(
+    code = "CLi = CL; Vi = V; dAdt[1] = -(CLi/Vi)*A[1]; CONC = A[1]/Vi",
+    state_init = "A[1] = TDM_INIT * Vi",
+    parameters = list(CL = 7.67, V = 97.7, TDM_INIT = 500),
+    obs = list(cmt = 1, scale = "Vi"),
+    declare_variables = c("CONC", "CLi", "Vi"),
+    cpp_show_code = FALSE
+  )
+
+  # --- From test_reparametrization.R: Carreno 2cmt model ---
+  covs_carreno <- list(
+    CRCL = new_covariate(5),
+    CL_HEMO = new_covariate(0)
+  )
+  model_carreno <- new_ode_model(
+    code = "
+    CLi = SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO \
+    Vi = V \
+    Qi = K12 * Vi \
+    V2i = Qi / K21 \
+    dAdt[0] = -(CLi/V)*A[0] - K12*A[0] + K21*A[1] \
+    dAdt[1] = K12*A[0] - K21*A[1] \
+    dAdt[2] = A[0]/V
+  ",
+    parameters = list(
+      V = 25.76, SCLSlope = 0.036, K12 = 2.29, K21 = 1.44,
+      SCLInter = 0.18, TDM_INIT = 0
+    ),
+    declare_variables = c("CLi", "Qi", "Vi", "V2i"),
+    covariates = covs_carreno,
+    obs = list(cmt = 1, scale = "V"),
+    reparametrization = list(
+      "CL" = "SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO",
+      "V" = "V",
+      "Q" = "K12 * V",
+      "V2" = "(K12 * V) / K21"
+    )
+  )
+}
diff --git a/tests/testthat/test_advan_with_auc.R b/tests/testthat/test_advan_with_auc.R
index ced7ef98..0cc9e8a5 100644
--- a/tests/testthat/test_advan_with_auc.R
+++ b/tests/testthat/test_advan_with_auc.R
@@ -1,43 +1,14 @@
+# Uses models and parameters defined in setup.R (conditional, NOT_CRAN only):
+# - mod_1cmt_auc, mod_2cmt_auc, mod_3cmt_auc
+# - parameters_advan_auc
+
 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,
@@ -47,7 +18,7 @@ if (identical(Sys.getenv("NOT_CRAN"), "true")) {
   )
   data_bolus <- advan_create_data(
     reg_bolus,
-    parameters = parameters,
+    parameters = parameters_advan_auc,
     cmts = 5,
     t_obs = t_obs
   )
@@ -61,7 +32,7 @@ if (identical(Sys.getenv("NOT_CRAN"), "true")) {
   )
   data_infusion <- advan_create_data(
     reg_infusion,
-    parameters = parameters,
+    parameters = parameters_advan_auc,
     cmts = 6,
     t_obs = t_obs
   )
@@ -71,7 +42,7 @@ 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)
+  res1_iv_ode <- sim(ode = mod_1cmt_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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))
@@ -84,7 +55,7 @@ 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)
+  res2_iv_ode <- sim(ode = mod_2cmt_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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)
@@ -106,7 +77,7 @@ 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)
+  res2_inf_ode <- sim(ode = mod_2cmt_auc, regimen = reg_infusion, parameters = parameters_advan_auc, t_obs = t_obs)
 
   expect_equal(
     round(res2_inf_r[res2_inf_r$TIME %in% t_obs,]$AUC, 5),
@@ -131,7 +102,7 @@ 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)
+  res3_iv_ode <- sim(ode = mod_3cmt_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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)
@@ -153,7 +124,7 @@ 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)
+  res3_iv_ode <- sim(ode = mod_3cmt_auc, regimen = reg_infusion, parameters = parameters_advan_auc, 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)
@@ -170,4 +141,3 @@ test_that("Three compartment iv ADVAN runs", {
     round(res3_iv_ode[res3_iv_ode$comp == 4,]$y, 5)
   )
 })
-
diff --git a/tests/testthat/test_calc_ss_analytic.R b/tests/testthat/test_calc_ss_analytic.R
index 43634725..54ca63c5 100644
--- a/tests/testthat/test_calc_ss_analytic.R
+++ b/tests/testthat/test_calc_ss_analytic.R
@@ -1,3 +1,7 @@
+# Uses models defined in setup.R:
+# - mod_1cmt_iv, mod_2cmt_iv, mod_1cmt_oral
+# - pk_3cmt_iv, pk_2cmt_oral, pk_3cmt_oral (conditional, NOT_CRAN only)
+
 # shared parameters
 dose <- 100
 interval <- 12
@@ -6,13 +10,6 @@ reg_oral <- new_regimen(amt = dose, interval = interval, n = n_ss, type = "oral"
 reg_bolus <- new_regimen(amt = dose, interval = interval, n = n_ss, type = "bolus")
 reg_inf <- new_regimen(amt = dose, interval = interval, n = n_ss, type = "infusion")
 t_obs <- max(reg_oral$dose_times) + interval
-# Uses models defined in setup.R
-
-if (identical(Sys.getenv("NOT_CRAN"), "true")) {
-  pk_3cmt_iv <- new_ode_model("pk_3cmt_iv")
-}
-
-#delta <- function(x, ref) { abs(x-ref)/ref }
 
 test_that("1 cmt oral", {
   par <- list(CL = 5, V = 100, KA = 1)
@@ -72,7 +69,6 @@ test_that("1-cmt iv infusion", {
 
 test_that("2-cmt oral", {
   skip_on_cran()
-  pk_2cmt_oral <- new_ode_model("pk_2cmt_oral")
   par <- list(CL = 5, V = 100, Q = 3, V2 = 150, KA = 1)
   res_ana <- calc_ss_analytic(f = "2cmt_oral", dose = dose, interval = interval, parameters = par)
   res_ode <- sim(pk_2cmt_oral, parameters = par, regimen = reg_oral, t_obs = t_obs, only_obs = F)$y
@@ -98,7 +94,6 @@ test_that("2-cmt infusion", {
 
 test_that("3-cmt oral", {
   skip_on_cran()
-  pk_3cmt_oral <- new_ode_model("pk_3cmt_oral")
   par <- list(CL = 5, V = 100, Q = 3, V2 = 150, Q2 = 6, V3 = 250, KA = 1)
   res_ana <- calc_ss_analytic(f = "3cmt_oral", dose = dose, interval = interval, parameters = par)
   res_ode <- sim(pk_3cmt_oral, parameters = par, regimen = reg_oral, t_obs = t_obs, only_obs = F)$y
diff --git a/tests/testthat/test_cmt_mapping.R b/tests/testthat/test_cmt_mapping.R
index 48af64e6..dfeea738 100644
--- a/tests/testthat/test_cmt_mapping.R
+++ b/tests/testthat/test_cmt_mapping.R
@@ -1,12 +1,9 @@
-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"),
-  cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
-)
+# Uses model defined in setup.R:
+# - pk1cmt_oral_cmt_mapping
 
 test_that("Compartment mapping is added to attributes", {
-  expect_equal(attr(pk1cmt_oral_code, "cmt_mapping")[["oral"]], 1)
-  expect_equal(attr(pk1cmt_oral_code, "cmt_mapping")[["infusion"]], 2)
+  expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["oral"]], 1)
+  expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["infusion"]], 2)
 })
 
 test_that("Admin route is interpreted and simulated correctly", {
@@ -18,7 +15,7 @@ test_that("Admin route is interpreted and simulated correctly", {
   )
   p <- list(KA = 1, CL = 5, V = 50)
   res <-  sim_ode(
-    ode = pk1cmt_oral_code,
+    ode = pk1cmt_oral_cmt_mapping,
     parameters = p,
     regimen = regimen,
     only_obs = FALSE
diff --git a/tests/testthat/test_compare_results.R b/tests/testthat/test_compare_results.R
index 73652e9f..c6746311 100644
--- a/tests/testthat/test_compare_results.R
+++ b/tests/testthat/test_compare_results.R
@@ -1,30 +1,17 @@
-## models: shared between tests and take a while to compile
-#  - oral models
-## Uses model defined in setup.R
+# Uses models defined in setup.R:
+# - mod_1cmt_oral
+# - mod_1cmt_iv
+# - dose_in_cmt_2
+
 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
+  obs = list(cmt = 2, scale = "V")
 )
 
-
 test_that("Library and custom C++ and code matches analytic soln", {
   p <- list(KA = 1, CL = 5, V = 50)
   t_obs <- c(0:72)
@@ -43,7 +30,7 @@ test_that("Library and custom C++ and code matches analytic soln", {
     only_obs=TRUE
   )
 
-  pk1cmt_oral_code <- sim_ode(
+  pk1cmt_oral_code_res <- sim_ode(
     ode = pk1cmt_oral_code,
     parameters = p,
     duplicate_t_obs = TRUE,
@@ -55,7 +42,7 @@ test_that("Library and custom C++ and code matches analytic soln", {
 
   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))
+  expect_equal(round(pk1cmt_oral_code_res$y, 3), round(pk1cmt_oral_anal_res, 3))
 })
 
 
@@ -87,6 +74,7 @@ test_that("test bug EmCo 20150925", {
 })
 
 test_that("model size is appropriate (bug: JeHi 20151204)", {
+  skip_on_cran()
   pk3cmt <- new_ode_model(
     code = "
       dAdt[1] = -KA*A[1];
diff --git a/tests/testthat/test_get_var_y.R b/tests/testthat/test_get_var_y.R
index 264843f6..fb130274 100644
--- a/tests/testthat/test_get_var_y.R
+++ b/tests/testthat/test_get_var_y.R
@@ -1,6 +1,9 @@
+# Uses models defined in setup.R:
+# - mod_1cmt_iv
+# - mod_2cmt_iv
+# - mod_1cmt_iv_mm (conditional, NOT_CRAN only)
 
 ## Set up simulations to test variance:
-## Uses model defined in setup.R
 reg <- new_regimen(
   amt = 100,
   n = 3,
@@ -130,13 +133,12 @@ test_that("Two compartment model", {
 
 test_that("One compartment with MM kinetics", {
   skip_on_cran()
-  mod3 <- new_ode_model("pk_1cmt_iv_mm")
   par3 <- list(VMAX = 5, KM = 5, V = 10)
   omega3 <- c(0.1,
               0.05, 0.1,
               0.01, 0.01, 0.1)
   res <- sim_ode(
-    mod3,
+    mod_1cmt_iv_mm,
     parameters = par3,
     t_obs = t_obs,
     regimen = reg,
@@ -144,14 +146,14 @@ test_that("One compartment with MM kinetics", {
   )
 
   v1 <- get_var_y(
-    model = mod3,
+    model = mod_1cmt_iv_mm,
     parameters = par3,
     t_obs = t_obs,
     regimen = reg,
     omega = omega3
   )
   v2 <- get_var_y(
-    model = mod3,
+    model = mod_1cmt_iv_mm,
     parameters = par3,
     t_obs = t_obs,
     regimen = reg,
diff --git a/tests/testthat/test_iov.R b/tests/testthat/test_iov.R
index 650d1388..e2241005 100644
--- a/tests/testthat/test_iov.R
+++ b/tests/testthat/test_iov.R
@@ -1,44 +1,10 @@
-pars <- list(
-  "kappa_CL_1" = 0,
-  "kappa_CL_2" = 0,
-  "kappa_CL_3" = 0,
-  "eta_CL" = 0,
-  "CL" = 5,
-  "V" = 50,
-  "KA" = 1
-)
-pars0 <- list(
-  "CL" = 5,
-  "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(
+# Uses models and parameters defined in setup.R:
+# - pk_iov_none (pk0)
+# - pk_iov (pk1)
+# - pars_iov
+# - pars_iov_no_iov
+
+reg_iov <- new_regimen(
   amt = 100,
   interval = 24,
   n = 5,
@@ -49,9 +15,9 @@ 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,
-      parameters = pars0,
-      regimen = reg1,
+      ode = pk_iov_none,
+      parameters = pars_iov_no_iov,
+      regimen = reg_iov,
       omega = c(
         0.3 # IIV in CL
       ),
@@ -67,9 +33,9 @@ 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,
-      parameters = pars,
-      regimen = reg1,
+      ode = pk_iov,
+      parameters = pars_iov,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
@@ -89,9 +55,9 @@ 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,
-      parameters = pars,
-      regimen = reg1,
+      ode = pk_iov,
+      parameters = pars_iov,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
@@ -113,9 +79,9 @@ test_that("IOV is added to parameters", {
   set.seed(32)
 
   dat <- sim(
-    ode = pk1,
-    parameters = pars,
-    regimen = reg1,
+    ode = pk_iov,
+    parameters = pars_iov,
+    regimen = reg_iov,
     omega = c(
       iov_var, # IOV in CL
       0, iov_var,
@@ -149,17 +115,18 @@ test_that("IOV is added to parameters", {
 })
 
 test_that("Change in F in 2nd bin is applied in 2nd bin and not later.", {
+  skip_on_cran()
   # Previously this was an issue because F, when defined in pk_code(), was not updated before the
   # dose was applied to the state vector, so the bioavailability was not applied at the right time.
   # This was fixed by rearranging the order of execution in sim.cpp in the main loop.
 
-  pars <- list(
+  pars_iov_f <- list(
     "CL" = 5,
     "V" = 50,
     "KA" = 1,
     "F" = 1
   )
-  pk1 <- new_ode_model(
+  pk_iov_f <- new_ode_model(
     code = "
         dAdt[1] = -KA * A[1]
         dAdt[2] = +KA * A[1] -(CLi/V) * A[2]
@@ -175,18 +142,18 @@ test_that("Change in F in 2nd bin is applied in 2nd bin and not later.", {
     obs = list(cmt = 2, scale = "V"),
     dose = list(cmt = 1, bioav = "Fi"),
     declare_variables = c("kappa_F", "Fi", "CLi"),
-    parameters = names(pars),
-    cpp_show_code = F
+    parameters = names(pars_iov_f),
+    cpp_show_code = FALSE
   )
   reg <- new_regimen(amt = 800, interval = 24, n = 10, type = "oral")
 
   # For a first simulation, we're simulating with no variability across the IOV bins:
-  pars$kappa_F_1 <- 0
-  pars$kappa_F_2 <- 0
-  pars$kappa_F_3 <- 0
+  pars_iov_f$kappa_F_1 <- 0
+  pars_iov_f$kappa_F_2 <- 0
+  pars_iov_f$kappa_F_3 <- 0
   args_sim1 <- args <- list(
-    ode = pk1,
-    parameters = pars,
+    ode = pk_iov_f,
+    parameters = pars_iov_f,
     regimen = reg,
     only_obs = TRUE,
     t_obs = seq(0, 50, .25),
@@ -194,10 +161,10 @@ test_that("Change in F in 2nd bin is applied in 2nd bin and not later.", {
   )
   # For a second simulation, we're applying a change in parameter for the 2nd bin (24-48 hrs).
   # This should affect predictions from 24 onward.
-  pars$kappa_F_2 <- 1 # 2nd bin
+  pars_iov_f$kappa_F_2 <- 1 # 2nd bin
   args_sim2 <- args <- list(
-    ode = pk1,
-    parameters = pars,
+    ode = pk_iov_f,
+    parameters = pars_iov_f,
     regimen = reg,
     only_obs = TRUE,
     t_obs = seq(0, 50, .25),
@@ -217,9 +184,9 @@ 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 = pk_iov,
       parameters_table = parameters_table,
-      regimen = reg1,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
@@ -238,10 +205,10 @@ test_that("error is not invoked when using parameters_table", {
   # specifying both parameters and parameters_table should fail
   expect_error(
     dat <- sim(
-      ode = pk1,
-      parameters = pars,
+      ode = pk_iov,
+      parameters = pars_iov,
       parameters_table = parameters_table,
-      regimen = reg1,
+      regimen = reg_iov,
       omega = c(
         iov_var, # IOV in CL
         0, iov_var,
diff --git a/tests/testthat/test_mixture_model.R b/tests/testthat/test_mixture_model.R
index 4289592a..447239b9 100644
--- a/tests/testthat/test_mixture_model.R
+++ b/tests/testthat/test_mixture_model.R
@@ -1,22 +1,15 @@
-## Setup model + params
-covs <- list(WT = PKPDsim::new_covariate(70))
-mod <- new_ode_model(
-  code = "
-      dAdt[0] = -(CL*(WT/70.0)/V)*A[0];
-    ",
-  pk_code = " ",
-  obs = list(cmt = 1, scale = "V"),
-  mixture = list(CL = list(values = c(5, 15), probability = 0.3)),
-  covariates = covs,
-)
-par <- list(CL = 3, V = 50)
-reg <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
-t_obs <- seq(0, 36, 4)
+# Uses model and covariates defined in setup.R:
+# - mod_mixture
+# - covs_mixture
+
+par_mixture <- list(CL = 3, V = 50)
+reg_mixture <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
+t_obs_mixture <- seq(0, 36, 4)
 
 test_that("mixture model works properly for single patient", {
-  res0 <- sim_ode(mod, parameters = par, regimen = reg, covariates = covs, t_obs = t_obs, only_obs=T) # mixture_group not supplied
-  res1 <- sim(mod, parameters = par, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = 1, only_obs=T)
-  res2 <- sim(mod, parameters = par, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = 2, only_obs=T)
+  res0 <- sim_ode(mod_mixture, parameters = par_mixture, regimen = reg_mixture, covariates = covs_mixture, t_obs = t_obs_mixture, only_obs=T) # mixture_group not supplied
+  res1 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T)
+  res2 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 2, only_obs=T)
   expect_equal(round(res0[res0$t == 24,]$y, 2), 9.07) # should use whatever is in `parameters`
   expect_equal(round(res1[res1$t == 24,]$y, 2), 5.82)
   expect_equal(round(res2[res2$t == 24,]$y, 2), 1.15)
@@ -25,9 +18,9 @@ test_that("mixture model works properly for single patient", {
 test_that("mixture model works properly when vectorized (using parameters_table)", {
   partab <- data.frame(CL = rep(0, 6), V = rep(50, 6))
   suppressMessages({
-    expect_error(sim_ode(mod, parameters_table = partab, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = 1, only_obs=T))
-    res1 <- sim(mod, parameters_table = partab, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = rep(1, 6), only_obs=T)
-    res2 <- sim(mod, parameters_table = partab, regimen = reg, t_obs = t_obs, covariates = covs, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
+    expect_error(sim_ode(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T))
+    res1 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(1, 6), only_obs=T)
+    res2 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
   })
   expect_equal(round(res1[res1$t == 24,]$y, 2), rep(5.82, 6))
   expect_equal(round(res2[res2$t == 24,]$y, 2), rep(c(5.82, 1.15), 3))
@@ -39,8 +32,8 @@ test_that("mixture model works properly when vectorized (using parameters_table)
 test_that("mixture model works properly when vectorized (using covariates_table)", {
   covtab <- data.frame(ID = 1:8, WT = rep(seq(40, 130, 30), 2))
   suppressMessages({
-    expect_error(sim(mod, parameters = par, covariates_table = covtab, regimen = reg, t_obs = t_obs, mixture_group = 1, only_obs=T))
-    res <- sim(mod, parameters = par, covariates_table = covtab, regimen = reg, t_obs = t_obs, mixture_group = rep(c(1, 2), each=4), only_obs=T)
+    expect_error(sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = 1, only_obs=T))
+    res <- sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = rep(c(1, 2), each=4), only_obs=T)
   })
   expect_equal(round(res[res$t == 24,]$y, 2), c(9.39, 5.82, 3.83, 2.65, 2.99, 1.15, 0.52, 0.25))
 })
diff --git a/tests/testthat/test_multi_obs.R b/tests/testthat/test_multi_obs.R
index f64511f7..15be1909 100644
--- a/tests/testthat/test_multi_obs.R
+++ b/tests/testthat/test_multi_obs.R
@@ -1,24 +1,17 @@
-# Example multiple observation types (e.g. different compartments! not just different residual errors)
+# Uses model defined in setup.R:
+# - pk_multi_obs
+# - vars_multi_obs
 
-## define parameters
-vars <- c("CONC", "METAB", "METAB2", "ACT")
-pk1 <- new_ode_model(
-  code = "dAdt[1] = -(CL/V)*A[1]; CONC = 1000*A[1]/V; METAB = CONC/2; METAB2 = CONC * t; ACT = 15",
-  obs = list(variable = vars),
-  declare_variables = vars,
-  cpp_show_code = F
-)
-
-regimen  <- new_regimen(amt = 100, interval = 12, n = 5, type="infusion", t_inf = 1)
-parameters   <- list("CL" = 15, "V" = 150)
-omega <- PKPDsim::cv_to_omega(list("CL" = 0.2, "V" = 0.2), parameters[1:2])
+regimen_multi_obs <- new_regimen(amt = 100, interval = 12, n = 5, type="infusion", t_inf = 1)
+parameters_multi_obs <- list("CL" = 15, "V" = 150)
+omega_multi_obs <- PKPDsim::cv_to_omega(list("CL" = 0.2, "V" = 0.2), parameters_multi_obs[1:2])
 
 test_that("obs types are output by `sim`", {
   obs_type <- c(1,2,1,3,1)
   data <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -33,9 +26,9 @@ test_that("obs types are output by `sim`", {
 test_that("check obs at same timepoint but with different obs_type", {
   obs_type <- c(1,2,1,3,1)
   t_same <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -58,9 +51,9 @@ test_that("check that residual error correctly applied to right var", {
   ruv_term1 <- list(prop = c(.1, 0, 0), add = c(1, 0, 0))
 
   data2 <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -75,9 +68,9 @@ test_that("check that residual error correctly applied to right var", {
   expect_true(data2$y[-4][4] != y[4])
 
   data3 <- sim(
-    ode = pk1,
+    ode = pk_multi_obs,
     parameters = list(CL = 20, V = 200),
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     int_step_size = 0.1,
     only_obs = TRUE,
     obs_type = obs_type,
@@ -99,10 +92,10 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
   ruv_multi <- list(prop = c(0.1, 1), add = c(1, 20))
 
   s_single <- sim(
-    ode = pk1,
-    parameters = parameters,
+    ode = pk_multi_obs,
+    parameters = parameters_multi_obs,
     n_ind = 1,
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     only_obs = TRUE,
     t_obs = c(2,4,6,8),
     seed = 123456,
@@ -112,10 +105,10 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
   ## specified as multi, but obs_type is all 1, so should
   ## give same results as s_single
   s_multi1 <- sim(
-    ode = pk1,
-    parameters = parameters,
+    ode = pk_multi_obs,
+    parameters = parameters_multi_obs,
     n_ind = 1,
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     only_obs = TRUE,
     t_obs = c(2,4,6,8),
     obs_type = c(1, 1, 1, 1),
@@ -128,10 +121,10 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
   ## specifying ruv as multi-type when multiple obs_type
   ## should now give different results
   s_multi1 <- sim(
-    ode = pk1,
-    parameters = parameters,
+    ode = pk_multi_obs,
+    parameters = parameters_multi_obs,
     n_ind = 1,
-    regimen = regimen,
+    regimen = regimen_multi_obs,
     only_obs = TRUE,
     t_obs = c(2,4,6,8),
     obs_type = c(1, 2, 1, 2),
@@ -145,9 +138,9 @@ test_that("specifying ruv as multi-type when only 1 obs_type", {
 
 test_that("multi-obs with baseline and obs_time = dose_time works correctly", {
   tmp <- sim(
-    pk1,
-    parameters = parameters,
-    regimen = regimen,
+    pk_multi_obs,
+    parameters = parameters_multi_obs,
+    regimen = regimen_multi_obs,
     t_obs = c(0, 0, 6, 6),
     obs_type = c(1, 4, 1, 4),
     only_obs = TRUE,
diff --git a/tests/testthat/test_reparametrization.R b/tests/testthat/test_reparametrization.R
index 747d3734..f3352def 100644
--- a/tests/testthat/test_reparametrization.R
+++ b/tests/testthat/test_reparametrization.R
@@ -1,3 +1,9 @@
+# Uses model and covariates defined in setup.R (conditional, NOT_CRAN only):
+# - model_carreno
+# - covs_carreno
+
+skip_on_cran()
+
 par_orig <- list(
   V = 25.76,
   SCLSlope = 0.036,
@@ -6,33 +12,8 @@ par_orig <- list(
   SCLInter = 0.18,
   TDM_INIT = 0)
 
-covs <- list(
-  CRCL = PKPDsim::new_covariate(5),
-  CL_HEMO = PKPDsim::new_covariate(0)
-)
-model <- new_ode_model( # Carreno et al
-  code = "
-    CLi = SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO \
-    Vi = V \
-    Qi = K12 * Vi \
-    V2i = Qi / K21 \
-    dAdt[0] = -(CLi/V)*A[0] - K12*A[0] + K21*A[1] \
-    dAdt[1] = K12*A[0] - K21*A[1] \
-    dAdt[2] = A[0]/V
-  ",
-  parameters = par_orig,
-  declare_variables = c("CLi", "Qi", "Vi", "V2i"),
-  covariates= covs,
-  obs = list(cmt = 1, scale = "V"),
-  reparametrization = list(
-    "CL" = "SCLInter + SCLSlope * (CRCL*16.667) + CL_HEMO",
-    "V" = "V",
-    "Q" = "K12 * V",
-    "V2" = "(K12 * V) / K21"
-  )
-)
-pars_covs_comb <- join_cov_and_par(covs, par_orig)
-pars <- reparametrize(model, pars_covs_comb, covariates = covs)
+pars_covs_comb <- join_cov_and_par(covs_carreno, par_orig)
+pars <- reparametrize(model_carreno, pars_covs_comb, covariates = covs_carreno)
 
 reg <- new_regimen(
   amt = 1000,
@@ -44,10 +25,10 @@ reg <- new_regimen(
   n_ss = 50
 )
 s <- sim(
-  ode = model,
+  ode = model_carreno,
   parameters = par_orig,
   regimen = reg,
-  covariates = covs,
+  covariates = covs_carreno,
   t_obs = c(0,24)
 )
 
@@ -96,4 +77,3 @@ test_that("Reparametrized model and analytics equations match", {
   cmin_ss_ode <- s[s$comp == "obs",]$y[1]
   expect_equal(cmin_ss_ode, cmin_ss_lin)
 })
-
diff --git a/tests/testthat/test_t_init.R b/tests/testthat/test_t_init.R
index f0e1a63a..abd4e2e5 100644
--- a/tests/testthat/test_t_init.R
+++ b/tests/testthat/test_t_init.R
@@ -1,23 +1,19 @@
-# Test t_init functionality
+# Uses model defined in setup.R (conditional, NOT_CRAN only):
+# - mod_t_init
+
+skip_on_cran()
 
+# Test t_init functionality
 ## e.g. TDM before first dose:
 ## at t=-8, conc=10000
 ## Use this as true value in the simulations
-par   <- list(CL = 7.67, V = 97.7, TDM_INIT = 500)
-mod <- new_ode_model(
-  code = "CLi = CL; Vi = V; dAdt[1] = -(CLi/Vi)*A[1]; CONC = A[1]/Vi",
-  state_init = "A[1] = TDM_INIT * Vi",
-  parameters = par,
-  obs = list(cmt = 1, scale = "Vi"),
-  declare_variables = c("CONC", "CLi", "Vi"),
-  cpp_show_code = F
-)
-reg <- new_regimen(amt = 100000, times=c(0, 24), type="bolus")
+par_t_init <- list(CL = 7.67, V = 97.7, TDM_INIT = 500)
+reg_t_init <- new_regimen(amt = 100000, times=c(0, 24), type="bolus")
 s <- sim(
-  ode = mod,
-  parameters = par,
+  ode = mod_t_init,
+  parameters = par_t_init,
   n_ind = 1,
-  regimen = reg,
+  regimen = reg_t_init,
   only_obs = TRUE,
   output_include = list(variables = TRUE),
   t_init = 10
@@ -33,5 +29,3 @@ test_that("TDM before first dose is considered a true initial value", {
 test_that("Variables are set (also in first row) when TDM before first dose", {
   expect_equal(round(s$CONC[1:5], 1), c(500, 462.2, 427.3, 395.1, 365.3))
 })
-
-
diff --git a/tests/testthat/test_timevar_cov.R b/tests/testthat/test_timevar_cov.R
index 60866dbc..b9b9248a 100644
--- a/tests/testthat/test_timevar_cov.R
+++ b/tests/testthat/test_timevar_cov.R
@@ -1,23 +1,16 @@
+# Uses model defined in setup.R (conditional, NOT_CRAN only):
+# - mod_2cmt_timevar
+
 skip_on_cran()
 
-par <- list(CL = 3, V = 50, Q = 2.5, V2 = 70)
-reg <- new_regimen(
+par_timevar <- list(CL = 3, V = 50, Q = 2.5, V2 = 70)
+reg_timevar <- new_regimen(
   amt = 250,
   n = 60,
   interval = 6,
   type = 'infusion',
   t_inf = 1
 )
-mod <- new_ode_model(
-  code = "
-      dAdt[1] = -(Q/V)*A[1] + (Q/V2)*A[2] -(CLi/V)*A[1];
-      dAdt[2] = -(Q/V2)*A[2] + (Q/V)*A[1];
-    ",
-  pk_code = "CLi = CL + CRCL",
-  obs = list(cmt = 2, scale = "V"),
-  covariates = list(CRCL = new_covariate(5)), declare_variables = "CLi",
-  cpp = FALSE
-)
 
 test_that("timevarying covariates handled", {
   # CLi changes by several orders of magnitude after
@@ -32,9 +25,9 @@ test_that("timevarying covariates handled", {
   )
   t_obs <- seq(0, 360, 0.1)
   sim1 <- sim_ode(
-    mod,
-    parameters = par,
-    regimen = reg,
+    mod_2cmt_timevar,
+    parameters = par_timevar,
+    regimen = reg_timevar,
     covariates = covs,
     only_obs = TRUE,
     t_obs = t_obs,
@@ -62,18 +55,18 @@ test_that("timevarying covariates are interpolated and affect PK", {
   )
   t_obs <- seq(0, 120, 2)
   sim2_inter <- sim_ode(
-    mod,
-    parameters = par,
-    regimen = reg,
+    mod_2cmt_timevar,
+    parameters = par_timevar,
+    regimen = reg_timevar,
     covariates = covs_inter,
     only_obs = TRUE,
     t_obs = t_obs,
     output_include = list(parameters = TRUE, covariates = TRUE, variables = TRUE)
   )
   sim2_locf <- sim_ode(
-    mod,
-    parameters = par,
-    regimen = reg,
+    mod_2cmt_timevar,
+    parameters = par_timevar,
+    regimen = reg_timevar,
     covariates = covs_locf,
     only_obs = TRUE,
     t_obs = t_obs,

From 9b1c4fbd1114781887af106df23348a5e1b15fc3 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 09:30:10 -0800
Subject: [PATCH 02/10] more updates

---
 tests/testthat/test_compare_results.R  |  2 +-
 tests/testthat/test_parameters_table.R | 10 +---------
 2 files changed, 2 insertions(+), 10 deletions(-)

diff --git a/tests/testthat/test_compare_results.R b/tests/testthat/test_compare_results.R
index c6746311..61120b02 100644
--- a/tests/testthat/test_compare_results.R
+++ b/tests/testthat/test_compare_results.R
@@ -12,7 +12,7 @@ pk1cmt_oral_code <- new_ode_model(
   obs = list(cmt = 2, scale = "V")
 )
 
-test_that("Library and custom C++ and code matches analytic soln", {
+test_that("Model from library and custom C++ and code matches analytic solution", {
   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
diff --git a/tests/testthat/test_parameters_table.R b/tests/testthat/test_parameters_table.R
index 1f9af030..e274d7c1 100644
--- a/tests/testthat/test_parameters_table.R
+++ b/tests/testthat/test_parameters_table.R
@@ -5,19 +5,11 @@ test_that("Simulating with table of parameters works", {
     V = rnorm(10, 5, 0.5)
   )
 
-  mod <- new_ode_model(
-    code = "
-      dAdt[1] = -(CL/V) * A[1];
-    ",
-    obs = list(cmt = 1, scale = "V"),
-    covariates = list(WT = new_covariate(70)),
-    cpp_show_code=FALSE
-  )
   par <- list(CL = 5, V = 50)
   reg <- new_regimen(amt = 2000, interval = 24, type = "infusion")
   covariates = list(WT = new_covariate(1))
   res <- sim_ode(
-    ode = mod,
+    ode = oral_1cmt_allometric,
     parameters_table = parameters_table,
     covariates = covariates,
     regimen = reg,

From 5be534040adaf80a4cb3c4b73059e76f6919d269 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 09:46:32 -0800
Subject: [PATCH 03/10] move iov

---
 tests/testthat/test_iov.R           | 226 ----------------------------
 tests/testthat/test_new_ode_model.R |  62 ++++++++
 tests/testthat/test_sim.R           | 163 ++++++++++++++++++++
 3 files changed, 225 insertions(+), 226 deletions(-)
 delete mode 100644 tests/testthat/test_iov.R

diff --git a/tests/testthat/test_iov.R b/tests/testthat/test_iov.R
deleted file mode 100644
index e2241005..00000000
--- a/tests/testthat/test_iov.R
+++ /dev/null
@@ -1,226 +0,0 @@
-# Uses models and parameters defined in setup.R:
-# - pk_iov_none (pk0)
-# - pk_iov (pk1)
-# - pars_iov
-# - pars_iov_no_iov
-
-reg_iov <- new_regimen(
-  amt = 100,
-  interval = 24,
-  n = 5,
-  type = "infusion"
-)
-iov_var <- 0.3 ^ 2 # 30% IOV
-
-test_that("Throws error when `iov_bins` supplied but not present in model", {
-  expect_error({
-    sim(
-      ode = pk_iov_none,
-      parameters = pars_iov_no_iov,
-      regimen = reg_iov,
-      omega = c(
-        0.3 # IIV in CL
-      ),
-      n = 1,
-      iov_bins = c(0, 24, 48, 72, 999), # !!
-      omega_type = "normal",
-      only_obs = TRUE,
-      output_include = list(parameters = TRUE, variables = TRUE)
-    )
-  }, "No IOV implemented for this model")
-})
-
-test_that("Throws error when number of `iov_bins` is higher than allowed for model", {
-  expect_error({
-    sim(
-      ode = pk_iov,
-      parameters = pars_iov,
-      regimen = reg_iov,
-      omega = c(
-        iov_var, # IOV in CL
-        0, iov_var,
-        0, 0, iov_var,
-        0, 0, 0, iov_var,
-        0, 0, 0, 0, 0.3 # IIV in CL
-      ),
-      n = 1,
-      iov_bins = c(0, 24, 48, 72, 999), # one bin too many
-      omega_type = "normal",
-      only_obs = TRUE,
-      output_include = list(parameters = TRUE, variables = TRUE)
-    )
-  }, "Number of allowed IOV bins for model is lower")
-})
-
-test_that("Throws warning when number of `iov_bins` is lower than allowed for model", {
-  expect_warning({
-    sim(
-      ode = pk_iov,
-      parameters = pars_iov,
-      regimen = reg_iov,
-      omega = c(
-        iov_var, # IOV in CL
-        0, iov_var,
-        0, 0, iov_var,
-        0, 0, 0, iov_var,
-        0, 0, 0, 0, 0.3 # IIV in CL
-      ),
-      n = 1,
-      iov_bins = c(0, 24, 999), # one bin too few
-      omega_type = "normal",
-      only_obs = TRUE,
-      output_include = list(parameters = TRUE, variables = TRUE)
-    )}, "Number of allowed IOV bins for model is higher"
-  )
-})
-
-test_that("IOV is added to parameters", {
-  skip_on_cran()
-  set.seed(32)
-
-  dat <- sim(
-    ode = pk_iov,
-    parameters = pars_iov,
-    regimen = reg_iov,
-    omega = c(
-      iov_var, # IOV in CL
-      0, iov_var,
-      0, 0, iov_var,
-      0, 0, 0, iov_var,
-      0, 0, 0, 0, 0.3 # IIV in CL
-    ),
-    n = 5,
-    iov_bins = c(0, 24, 48, 72),
-    omega_type = "normal",
-    only_obs = TRUE,
-    output_include = list(parameters = TRUE, variables = TRUE)
-  )
-
-  expect_equal(
-    signif(dat$kappa_CL[dat$t == 12], 4),
-    signif(dat$kappa_CL_1[dat$t == 1], 4)
-  )
-  expect_equal(
-    signif(dat$kappa_CL[dat$t == 36], 4),
-    signif(dat$kappa_CL_2[dat$t == 1], 4)
-  )
-  expect_equal(
-    signif(dat$kappa_CL[dat$t == 48], 4),
-    signif(dat$kappa_CL_3[dat$t == 1], 4)
-  )
-  expect_equal(
-    signif(exp(dat$kappa_CL + dat$eta_CL) * dat$CL, 4),
-    signif(dat$CL_iov, 4)
-  )
-})
-
-test_that("Change in F in 2nd bin is applied in 2nd bin and not later.", {
-  skip_on_cran()
-  # Previously this was an issue because F, when defined in pk_code(), was not updated before the
-  # dose was applied to the state vector, so the bioavailability was not applied at the right time.
-  # This was fixed by rearranging the order of execution in sim.cpp in the main loop.
-
-  pars_iov_f <- list(
-    "CL" = 5,
-    "V" = 50,
-    "KA" = 1,
-    "F" = 1
-  )
-  pk_iov_f <- new_ode_model(
-    code = "
-        dAdt[1] = -KA * A[1]
-        dAdt[2] = +KA * A[1] -(CLi/V) * A[2]
-    ",
-    pk_code = "
-        Fi = F * exp(kappa_F);
-        CLi = CL;
-    ",
-    iov = list(
-      cv = list(F = 0.2),
-      n_bins = 3
-    ),
-    obs = list(cmt = 2, scale = "V"),
-    dose = list(cmt = 1, bioav = "Fi"),
-    declare_variables = c("kappa_F", "Fi", "CLi"),
-    parameters = names(pars_iov_f),
-    cpp_show_code = FALSE
-  )
-  reg <- new_regimen(amt = 800, interval = 24, n = 10, type = "oral")
-
-  # For a first simulation, we're simulating with no variability across the IOV bins:
-  pars_iov_f$kappa_F_1 <- 0
-  pars_iov_f$kappa_F_2 <- 0
-  pars_iov_f$kappa_F_3 <- 0
-  args_sim1 <- args <- list(
-    ode = pk_iov_f,
-    parameters = pars_iov_f,
-    regimen = reg,
-    only_obs = TRUE,
-    t_obs = seq(0, 50, .25),
-    iov_bins = c(0L, 24L, 48L, 9999L)
-  )
-  # For a second simulation, we're applying a change in parameter for the 2nd bin (24-48 hrs).
-  # This should affect predictions from 24 onward.
-  pars_iov_f$kappa_F_2 <- 1 # 2nd bin
-  args_sim2 <- args <- list(
-    ode = pk_iov_f,
-    parameters = pars_iov_f,
-    regimen = reg,
-    only_obs = TRUE,
-    t_obs = seq(0, 50, .25),
-    iov_bins = c(0L, 24L, 48L, 9999L)
-  )
-  res1 <- do.call("sim_ode", args = args_sim1)
-  res2 <- do.call("sim_ode", args = args_sim2)
-  expect_true(min(res1[res1$y != res2$y,]$t) <= 25)
-})
-
-test_that("error is not invoked when using parameters_table", {
-  parameters_table <- data.frame(
-    CL = runif(10), V = runif(10), KA = runif(10), eta_CL = runif(10),
-    kappa_CL_1 = 0, kappa_CL_2 = 0, kappa_CL_3 = 0, kappa_CL_4 = 0
-  )
-
-  # specifying both parameters_table but for a model with IOV should not fail!
-  expect_silent(
-    dat <- sim(
-      ode = pk_iov,
-      parameters_table = parameters_table,
-      regimen = reg_iov,
-      omega = c(
-        iov_var, # IOV in CL
-        0, iov_var,
-        0, 0, iov_var,
-        0, 0, 0, iov_var,
-        0, 0, 0, 0, 0.3 # IIV in CL
-      ),
-      n = 5,
-      iov_bins = c(0, 24, 48, 72),
-      omega_type = "normal",
-      only_obs = TRUE,
-      output_include = list(parameters = TRUE, variables = TRUE)
-    )
-  )
-
-  # specifying both parameters and parameters_table should fail
-  expect_error(
-    dat <- sim(
-      ode = pk_iov,
-      parameters = pars_iov,
-      parameters_table = parameters_table,
-      regimen = reg_iov,
-      omega = c(
-        iov_var, # IOV in CL
-        0, iov_var,
-        0, 0, iov_var,
-        0, 0, 0, iov_var,
-        0, 0, 0, 0, 0.3 # IIV in CL
-      ),
-      n = 5,
-      iov_bins = c(0, 24, 48, 72),
-      omega_type = "normal",
-      only_obs = TRUE,
-      output_include = list(parameters = TRUE, variables = TRUE)
-    )
-  )
-})
diff --git a/tests/testthat/test_new_ode_model.R b/tests/testthat/test_new_ode_model.R
index 5b5141a6..c27f3205 100644
--- a/tests/testthat/test_new_ode_model.R
+++ b/tests/testthat/test_new_ode_model.R
@@ -163,3 +163,65 @@ test_that("specifying overlapping covariates and variables throws error", {
     )
   })
 })
+
+test_that("IOV: Change in F in 2nd bin is applied in 2nd bin and not later.", {
+  skip_on_cran()
+  # Previously this was an issue because F, when defined in pk_code(), was not updated before the
+  # dose was applied to the state vector, so the bioavailability was not applied at the right time.
+  # This was fixed by rearranging the order of execution in sim.cpp in the main loop.
+
+  pars_iov_f <- list(
+    "CL" = 5,
+    "V" = 50,
+    "KA" = 1,
+    "F" = 1
+  )
+  pk_iov_f <- new_ode_model(
+    code = "
+        dAdt[1] = -KA * A[1]
+        dAdt[2] = +KA * A[1] -(CLi/V) * A[2]
+    ",
+    pk_code = "
+        Fi = F * exp(kappa_F);
+        CLi = CL;
+    ",
+    iov = list(
+      cv = list(F = 0.2),
+      n_bins = 3
+    ),
+    obs = list(cmt = 2, scale = "V"),
+    dose = list(cmt = 1, bioav = "Fi"),
+    declare_variables = c("kappa_F", "Fi", "CLi"),
+    parameters = names(pars_iov_f),
+    cpp_show_code = FALSE
+  )
+  reg <- new_regimen(amt = 800, interval = 24, n = 10, type = "oral")
+
+  # For a first simulation, we're simulating with no variability across the IOV bins:
+  pars_iov_f$kappa_F_1 <- 0
+  pars_iov_f$kappa_F_2 <- 0
+  pars_iov_f$kappa_F_3 <- 0
+  args_sim1 <- args <- list(
+    ode = pk_iov_f,
+    parameters = pars_iov_f,
+    regimen = reg,
+    only_obs = TRUE,
+    t_obs = seq(0, 50, .25),
+    iov_bins = c(0L, 24L, 48L, 9999L)
+  )
+  # For a second simulation, we're applying a change in parameter for the 2nd bin (24-48 hrs).
+  # This should affect predictions from 24 onward.
+  pars_iov_f$kappa_F_2 <- 1 # 2nd bin
+  args_sim2 <- args <- list(
+    ode = pk_iov_f,
+    parameters = pars_iov_f,
+    regimen = reg,
+    only_obs = TRUE,
+    t_obs = seq(0, 50, .25),
+    iov_bins = c(0L, 24L, 48L, 9999L)
+  )
+  res1 <- do.call("sim_ode", args = args_sim1)
+  res2 <- do.call("sim_ode", args = args_sim2)
+  expect_true(min(res1[res1$y != res2$y,]$t) <= 25)
+})
+
diff --git a/tests/testthat/test_sim.R b/tests/testthat/test_sim.R
index cb676754..6f80929c 100644
--- a/tests/testthat/test_sim.R
+++ b/tests/testthat/test_sim.R
@@ -578,3 +578,166 @@ test_that("t_max is shifted correctly when t_ss != 0", {
   ))
   expect_equal(sum(is.na(evtab3)), 0)
 })
+
+describe("IOV", {
+
+  reg_iov <- new_regimen(
+    amt = 100,
+    interval = 24,
+    n = 5,
+    type = "infusion"
+  )
+  iov_var <- 0.3 ^ 2 # 30% IOV
+
+  test_that("Throws error when `iov_bins` supplied but not present in model", {
+    expect_error({
+      sim(
+        ode = pk_iov_none,
+        parameters = pars_iov_no_iov,
+        regimen = reg_iov,
+        omega = c(
+          0.3 # IIV in CL
+        ),
+        n = 1,
+        iov_bins = c(0, 24, 48, 72, 999), # !!
+        omega_type = "normal",
+        only_obs = TRUE,
+        output_include = list(parameters = TRUE, variables = TRUE)
+      )
+    }, "No IOV implemented for this model")
+  })
+
+  test_that("Throws error when number of `iov_bins` is higher than allowed for model", {
+    expect_error({
+      sim(
+        ode = pk_iov,
+        parameters = pars_iov,
+        regimen = reg_iov,
+        omega = c(
+          iov_var, # IOV in CL
+          0, iov_var,
+          0, 0, iov_var,
+          0, 0, 0, iov_var,
+          0, 0, 0, 0, 0.3 # IIV in CL
+        ),
+        n = 1,
+        iov_bins = c(0, 24, 48, 72, 999), # one bin too many
+        omega_type = "normal",
+        only_obs = TRUE,
+        output_include = list(parameters = TRUE, variables = TRUE)
+      )
+    }, "Number of allowed IOV bins for model is lower")
+  })
+
+  test_that("Throws warning when number of `iov_bins` is lower than allowed for model", {
+    expect_warning({
+      sim(
+        ode = pk_iov,
+        parameters = pars_iov,
+        regimen = reg_iov,
+        omega = c(
+          iov_var, # IOV in CL
+          0, iov_var,
+          0, 0, iov_var,
+          0, 0, 0, iov_var,
+          0, 0, 0, 0, 0.3 # IIV in CL
+        ),
+        n = 1,
+        iov_bins = c(0, 24, 999), # one bin too few
+        omega_type = "normal",
+        only_obs = TRUE,
+        output_include = list(parameters = TRUE, variables = TRUE)
+      )}, "Number of allowed IOV bins for model is higher"
+    )
+  })
+
+  test_that("IOV is added to parameters", {
+    skip_on_cran()
+    set.seed(32)
+
+    dat <- sim(
+      ode = pk_iov,
+      parameters = pars_iov,
+      regimen = reg_iov,
+      omega = c(
+        iov_var, # IOV in CL
+        0, iov_var,
+        0, 0, iov_var,
+        0, 0, 0, iov_var,
+        0, 0, 0, 0, 0.3 # IIV in CL
+      ),
+      n = 5,
+      iov_bins = c(0, 24, 48, 72),
+      omega_type = "normal",
+      only_obs = TRUE,
+      output_include = list(parameters = TRUE, variables = TRUE)
+    )
+
+    expect_equal(
+      signif(dat$kappa_CL[dat$t == 12], 4),
+      signif(dat$kappa_CL_1[dat$t == 1], 4)
+    )
+    expect_equal(
+      signif(dat$kappa_CL[dat$t == 36], 4),
+      signif(dat$kappa_CL_2[dat$t == 1], 4)
+    )
+    expect_equal(
+      signif(dat$kappa_CL[dat$t == 48], 4),
+      signif(dat$kappa_CL_3[dat$t == 1], 4)
+    )
+    expect_equal(
+      signif(exp(dat$kappa_CL + dat$eta_CL) * dat$CL, 4),
+      signif(dat$CL_iov, 4)
+    )
+  })
+  
+  test_that("error is not invoked when using parameters_table", {
+    parameters_table <- data.frame(
+      CL = runif(10), V = runif(10), KA = runif(10), eta_CL = runif(10),
+      kappa_CL_1 = 0, kappa_CL_2 = 0, kappa_CL_3 = 0, kappa_CL_4 = 0
+    )
+
+    # specifying both parameters_table but for a model with IOV should not fail!
+    expect_silent(
+      dat <- sim(
+        ode = pk_iov,
+        parameters_table = parameters_table,
+        regimen = reg_iov,
+        omega = c(
+          iov_var, # IOV in CL
+          0, iov_var,
+          0, 0, iov_var,
+          0, 0, 0, iov_var,
+          0, 0, 0, 0, 0.3 # IIV in CL
+        ),
+        n = 5,
+        iov_bins = c(0, 24, 48, 72),
+        omega_type = "normal",
+        only_obs = TRUE,
+        output_include = list(parameters = TRUE, variables = TRUE)
+      )
+    )
+
+    # specifying both parameters and parameters_table should fail
+    expect_error(
+      dat <- sim(
+        ode = pk_iov,
+        parameters = pars_iov,
+        parameters_table = parameters_table,
+        regimen = reg_iov,
+        omega = c(
+          iov_var, # IOV in CL
+          0, iov_var,
+          0, 0, iov_var,
+          0, 0, 0, iov_var,
+          0, 0, 0, 0, 0.3 # IIV in CL
+        ),
+        n = 5,
+        iov_bins = c(0, 24, 48, 72),
+        omega_type = "normal",
+        only_obs = TRUE,
+        output_include = list(parameters = TRUE, variables = TRUE)
+      )
+    )
+  })
+})

From 252e27728d56840d12006afe1b4986bd4e1bd79c Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 09:52:14 -0800
Subject: [PATCH 04/10] consolidate more tests

---
 R/sim.R                                   | 100 +++++++++
 tests/testthat/test_check_mixture_model.R |  17 ++
 tests/testthat/test_cmt_mapping.R         |  53 -----
 tests/testthat/test_compare_results.R     | 254 ---------------------
 tests/testthat/test_mixture_model.R       |  57 -----
 tests/testthat/test_sim.R                 | 259 ++++++++++++++++++++++
 6 files changed, 376 insertions(+), 364 deletions(-)
 create mode 100644 tests/testthat/test_check_mixture_model.R
 delete mode 100644 tests/testthat/test_cmt_mapping.R
 delete mode 100644 tests/testthat/test_compare_results.R
 delete mode 100644 tests/testthat/test_mixture_model.R

diff --git a/R/sim.R b/R/sim.R
index feda5c67..b8eb9a26 100644
--- a/R/sim.R
+++ b/R/sim.R
@@ -653,3 +653,103 @@ sim <- function (ode = NULL,
   attr(comb, "parameters") <- attr(ode, "parameters")
   return(comb)
 }
+
+describe("Mixture models", {
+  # Uses model and covariates defined in setup.R:
+  # - mod_mixture
+  # - covs_mixture
+
+  par_mixture <- list(CL = 3, V = 50)
+  reg_mixture <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
+  t_obs_mixture <- seq(0, 36, 4)
+
+  test_that("mixture model works properly for single patient", {
+    res0 <- sim_ode(mod_mixture, parameters = par_mixture, regimen = reg_mixture, covariates = covs_mixture, t_obs = t_obs_mixture, only_obs=T) # mixture_group not supplied
+    res1 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T)
+    res2 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 2, only_obs=T)
+    expect_equal(round(res0[res0$t == 24,]$y, 2), 9.07) # should use whatever is in `parameters`
+    expect_equal(round(res1[res1$t == 24,]$y, 2), 5.82)
+    expect_equal(round(res2[res2$t == 24,]$y, 2), 1.15)
+  })
+
+  test_that("mixture model works properly when vectorized (using parameters_table)", {
+    partab <- data.frame(CL = rep(0, 6), V = rep(50, 6))
+    suppressMessages({
+      expect_error(sim_ode(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T))
+      res1 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(1, 6), only_obs=T)
+      res2 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
+    })
+    expect_equal(round(res1[res1$t == 24,]$y, 2), rep(5.82, 6))
+    expect_equal(round(res2[res2$t == 24,]$y, 2), rep(c(5.82, 1.15), 3))
+    expect_equal(res2[res2$id == 1,]$CL[1], 5)
+    expect_equal(res2[res2$id == 2,]$CL[1], 15)
+    expect_equal(res2[res2$id == 3,]$CL[1], 5)
+  })
+
+  test_that("mixture model works properly when vectorized (using covariates_table)", {
+    covtab <- data.frame(ID = 1:8, WT = rep(seq(40, 130, 30), 2))
+    suppressMessages({
+      expect_error(sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = 1, only_obs=T))
+      res <- sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = rep(c(1, 2), each=4), only_obs=T)
+    })
+    expect_equal(round(res[res$t == 24,]$y, 2), c(9.39, 5.82, 3.83, 2.65, 2.99, 1.15, 0.52, 0.25))
+  })
+
+})
+
+describe("Compartment mapping", {
+  # Uses model defined in setup.R:
+  # - pk1cmt_oral_cmt_mapping
+
+  test_that("Compartment mapping is added to attributes", {
+    expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["oral"]], 1)
+    expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["infusion"]], 2)
+  })
+
+  test_that("Admin route is interpreted and simulated correctly", {
+    regimen <- new_regimen(
+      amt = c(100, 100, 100, 100),
+      times = c(0, 12, 24, 36),
+      type = c("oral", "oral", "infusion", "infusion"),
+      t_inf = c(0, 0, 1, 1)
+    )
+    p <- list(KA = 1, CL = 5, V = 50)
+    res <-  sim_ode(
+      ode = pk1cmt_oral_cmt_mapping,
+      parameters = p,
+      regimen = regimen,
+      only_obs = FALSE
+    )
+    expect_equal(round(res$y[res$comp == 1 & res$t == 25], 4), 2e-04)
+    expect_true(res$y[res$comp == 2 & res$t == 25] >= 100)
+  })
+
+  test_that("multiple scaling types on one compartment works", {
+    skip_on_cran()
+    mod <- new_ode_model(
+      code = "
+        dAdt[1] = -KA * A[1];
+        dAdt[2] = -(CL/V) * A[2] + KA*A[1];
+      ",
+      obs = list(
+        cmt = c(2, 2),
+        scale = c(1, "V"),
+        label = c("abs", "conc")
+      ),
+      cpp_show_code = FALSE
+    )
+    par <- list(CL = 5, V = 50, KA = .5)
+    reg <- new_regimen(amt = 100, n = 5, interval = 12)
+    res <- sim_ode(
+      ode = mod,
+      parameters = par,
+      regimen = reg,
+      only_obs = TRUE
+    )
+    dat <- cbind(res[res$comp == "abs",]$y, res[res$comp == "conc",]$y)
+    expect_true("PKPDsim" %in% class(mod))
+    expect_equal(length(unique(res$comp)), 2)
+    expect_equal(round(dat[,1],1), round(dat[,2]*par$V,1))
+  })
+
+})
\ No newline at end of file
diff --git a/tests/testthat/test_check_mixture_model.R b/tests/testthat/test_check_mixture_model.R
new file mode 100644
index 00000000..8e97be4f
--- /dev/null
+++ b/tests/testthat/test_check_mixture_model.R
@@ -0,0 +1,17 @@
+test_that("expected mixture model errors are caught", {
+  mixture1 <- list(CL = list(values = c(5, 10), probability = 0.3))
+  mixture2 <- list(
+    CL = list(values = c(5, 10), probability = 0.3),
+    V = list(values = c(10, 20), probability = 0.5)
+  )
+  mixture_multi <- list(
+    CL = list(values = c(5, 6, 7), probability =c(0.3, 0.3, 0.4))
+  )
+  mixture_oops1 <- list(CL = list(values = c(5, 10), probability = 30))
+  mixture_oops2 <- list(CL = list(values = c(5, 10), probability = -0.3))
+  expect_error(check_mixture_model(mixture2, c("CL", "V")))
+  expect_error(check_mixture_model(mixture1, c("Ke", "V")))
+  expect_error(check_mixture_model(mixture_multi, c("CL", "V")))
+  expect_error(check_mixture_model(mixture_oops1, c("CL", "V")))
+  expect_error(check_mixture_model(mixture_oops2, c("CL", "V")))
+})
diff --git a/tests/testthat/test_cmt_mapping.R b/tests/testthat/test_cmt_mapping.R
deleted file mode 100644
index dfeea738..00000000
--- a/tests/testthat/test_cmt_mapping.R
+++ /dev/null
@@ -1,53 +0,0 @@
-# Uses model defined in setup.R:
-# - pk1cmt_oral_cmt_mapping
-
-test_that("Compartment mapping is added to attributes", {
-  expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["oral"]], 1)
-  expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["infusion"]], 2)
-})
-
-test_that("Admin route is interpreted and simulated correctly", {
-  regimen <- new_regimen(
-    amt = c(100, 100, 100, 100),
-    times = c(0, 12, 24, 36),
-    type = c("oral", "oral", "infusion", "infusion"),
-    t_inf = c(0, 0, 1, 1)
-  )
-  p <- list(KA = 1, CL = 5, V = 50)
-  res <-  sim_ode(
-    ode = pk1cmt_oral_cmt_mapping,
-    parameters = p,
-    regimen = regimen,
-    only_obs = FALSE
-  )
-  expect_equal(round(res$y[res$comp == 1 & res$t == 25], 4), 2e-04)
-  expect_true(res$y[res$comp == 2 & res$t == 25] >= 100)
-})
-
-test_that("multiple scaling types on one compartment works", {
-  skip_on_cran()
-  mod <- new_ode_model(
-    code = "
-      dAdt[1] = -KA * A[1];
-      dAdt[2] = -(CL/V) * A[2] + KA*A[1];
-    ",
-    obs = list(
-      cmt = c(2, 2),
-      scale = c(1, "V"),
-      label = c("abs", "conc")
-    ),
-    cpp_show_code = FALSE
-  )
-  par <- list(CL = 5, V = 50, KA = .5)
-  reg <- new_regimen(amt = 100, n = 5, interval = 12)
-  res <- sim_ode(
-    ode = mod,
-    parameters = par,
-    regimen = reg,
-    only_obs = TRUE
-  )
-  dat <- cbind(res[res$comp == "abs",]$y, res[res$comp == "conc",]$y)
-  expect_true("PKPDsim" %in% class(mod))
-  expect_equal(length(unique(res$comp)), 2)
-  expect_equal(round(dat[,1],1), round(dat[,2]*par$V,1))
-})
diff --git a/tests/testthat/test_compare_results.R b/tests/testthat/test_compare_results.R
deleted file mode 100644
index 61120b02..00000000
--- a/tests/testthat/test_compare_results.R
+++ /dev/null
@@ -1,254 +0,0 @@
-# Uses models defined in setup.R:
-# - mod_1cmt_oral
-# - mod_1cmt_iv
-# - dose_in_cmt_2
-
-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")
-)
-
-test_that("Model from library and custom C++ and code matches analytic solution", {
-  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 = 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_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)", {
-  skip_on_cran()
-  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_mixture_model.R b/tests/testthat/test_mixture_model.R
deleted file mode 100644
index 447239b9..00000000
--- a/tests/testthat/test_mixture_model.R
+++ /dev/null
@@ -1,57 +0,0 @@
-# Uses model and covariates defined in setup.R:
-# - mod_mixture
-# - covs_mixture
-
-par_mixture <- list(CL = 3, V = 50)
-reg_mixture <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
-t_obs_mixture <- seq(0, 36, 4)
-
-test_that("mixture model works properly for single patient", {
-  res0 <- sim_ode(mod_mixture, parameters = par_mixture, regimen = reg_mixture, covariates = covs_mixture, t_obs = t_obs_mixture, only_obs=T) # mixture_group not supplied
-  res1 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T)
-  res2 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 2, only_obs=T)
-  expect_equal(round(res0[res0$t == 24,]$y, 2), 9.07) # should use whatever is in `parameters`
-  expect_equal(round(res1[res1$t == 24,]$y, 2), 5.82)
-  expect_equal(round(res2[res2$t == 24,]$y, 2), 1.15)
-})
-
-test_that("mixture model works properly when vectorized (using parameters_table)", {
-  partab <- data.frame(CL = rep(0, 6), V = rep(50, 6))
-  suppressMessages({
-    expect_error(sim_ode(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T))
-    res1 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(1, 6), only_obs=T)
-    res2 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
-  })
-  expect_equal(round(res1[res1$t == 24,]$y, 2), rep(5.82, 6))
-  expect_equal(round(res2[res2$t == 24,]$y, 2), rep(c(5.82, 1.15), 3))
-  expect_equal(res2[res2$id == 1,]$CL[1], 5)
-  expect_equal(res2[res2$id == 2,]$CL[1], 15)
-  expect_equal(res2[res2$id == 3,]$CL[1], 5)
-})
-
-test_that("mixture model works properly when vectorized (using covariates_table)", {
-  covtab <- data.frame(ID = 1:8, WT = rep(seq(40, 130, 30), 2))
-  suppressMessages({
-    expect_error(sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = 1, only_obs=T))
-    res <- sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = rep(c(1, 2), each=4), only_obs=T)
-  })
-  expect_equal(round(res[res$t == 24,]$y, 2), c(9.39, 5.82, 3.83, 2.65, 2.99, 1.15, 0.52, 0.25))
-})
-
-test_that("expected mixture model errors are caught", {
-  mixture1 <- list(CL = list(values = c(5, 10), probability = 0.3))
-  mixture2 <- list(
-    CL = list(values = c(5, 10), probability = 0.3),
-    V = list(values = c(10, 20), probability = 0.5)
-  )
-  mixture_multi <- list(
-    CL = list(values = c(5, 6, 7), probability =c(0.3, 0.3, 0.4))
-  )
-  mixture_oops1 <- list(CL = list(values = c(5, 10), probability = 30))
-  mixture_oops2 <- list(CL = list(values = c(5, 10), probability = -0.3))
-  expect_error(check_mixture_model(mixture2, c("CL", "V")))
-  expect_error(check_mixture_model(mixture1, c("Ke", "V")))
-  expect_error(check_mixture_model(mixture_multi, c("CL", "V")))
-  expect_error(check_mixture_model(mixture_oops1, c("CL", "V")))
-  expect_error(check_mixture_model(mixture_oops2, c("CL", "V")))
-})
diff --git a/tests/testthat/test_sim.R b/tests/testthat/test_sim.R
index 6f80929c..9abb3789 100644
--- a/tests/testthat/test_sim.R
+++ b/tests/testthat/test_sim.R
@@ -741,3 +741,262 @@ describe("IOV", {
     )
   })
 })
+
+describe("Compare results from sims with references", {
+  
+  # Uses models defined in setup.R:
+  # - mod_1cmt_oral
+  # - mod_1cmt_iv
+  # - dose_in_cmt_2
+
+  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")
+  )
+
+  test_that("Model from library and custom C++ and code matches analytic solution", {
+    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 = 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_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)", {
+    skip_on_cran()
+    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)
+  })
+
+})
\ No newline at end of file

From 1250f664e2dce9ee9827c8b01b67e81a9f93e976 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 11:18:16 -0800
Subject: [PATCH 05/10] move tests

---
 R/{adherence.R => new_adherence.R}            |   0
 tests/testthat/test_advan.R                   | 210 ++++++++++++++++++
 tests/testthat/test_advan_with_auc.R          | 143 ------------
 tests/testthat/test_advan_with_covariates.R   |  62 ------
 ...{test_adherence.R => test_new_adherence.R} |   0
 5 files changed, 210 insertions(+), 205 deletions(-)
 rename R/{adherence.R => new_adherence.R} (100%)
 delete mode 100644 tests/testthat/test_advan_with_auc.R
 delete mode 100644 tests/testthat/test_advan_with_covariates.R
 rename tests/testthat/{test_adherence.R => test_new_adherence.R} (100%)

diff --git a/R/adherence.R b/R/new_adherence.R
similarity index 100%
rename from R/adherence.R
rename to R/new_adherence.R
diff --git a/tests/testthat/test_advan.R b/tests/testthat/test_advan.R
index 5a8dce3a..86435c06 100644
--- a/tests/testthat/test_advan.R
+++ b/tests/testthat/test_advan.R
@@ -124,3 +124,213 @@ test_that("Three compartment IV oral", {
   expect_true(!any(is.na(res3_oral$DV)))
   expect_equal(res3_oral, res3_oral_c)
 })
+
+describe("ADVANs with AUC", {
+  # Uses models and parameters defined in setup.R (conditional, NOT_CRAN only):
+  # - mod_1cmt_auc, mod_2cmt_auc, mod_3cmt_auc
+  # - parameters_advan_auc
+
+  if (identical(Sys.getenv("NOT_CRAN"), "true")) {
+    dose <- 100
+    interval <- 12
+    n_days <- 5
+    t_inf <- 1.5
+    t_obs <- c(3, 6, 8, 23, 47)
+
+    ## 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_advan_auc,
+      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_advan_auc,
+      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_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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_auc, regimen = reg_infusion, parameters = parameters_advan_auc, 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_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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_auc, regimen = reg_infusion, parameters = parameters_advan_auc, 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)
+    )
+  })
+
+})
+
+describe("ADVANs with covariates", {
+  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))
+  })
+})
\ No newline at end of file
diff --git a/tests/testthat/test_advan_with_auc.R b/tests/testthat/test_advan_with_auc.R
deleted file mode 100644
index 0cc9e8a5..00000000
--- a/tests/testthat/test_advan_with_auc.R
+++ /dev/null
@@ -1,143 +0,0 @@
-# Uses models and parameters defined in setup.R (conditional, NOT_CRAN only):
-# - mod_1cmt_auc, mod_2cmt_auc, mod_3cmt_auc
-# - parameters_advan_auc
-
-if (identical(Sys.getenv("NOT_CRAN"), "true")) {
-  dose <- 100
-  interval <- 12
-  n_days <- 5
-  t_inf <- 1.5
-  t_obs <- c(3, 6, 8, 23, 47)
-
-  ## 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_advan_auc,
-    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_advan_auc,
-    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_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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_auc, regimen = reg_infusion, parameters = parameters_advan_auc, 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_auc, regimen = reg_bolus, parameters = parameters_advan_auc, 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_auc, regimen = reg_infusion, parameters = parameters_advan_auc, 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_adherence.R b/tests/testthat/test_new_adherence.R
similarity index 100%
rename from tests/testthat/test_adherence.R
rename to tests/testthat/test_new_adherence.R

From be54964ad479d0e2fd5220f773f1a74baa7b4e75 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 11:20:47 -0800
Subject: [PATCH 06/10] fix test

---
 tests/testthat/test_parameters_table.R | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/testthat/test_parameters_table.R b/tests/testthat/test_parameters_table.R
index e274d7c1..4dcacdd1 100644
--- a/tests/testthat/test_parameters_table.R
+++ b/tests/testthat/test_parameters_table.R
@@ -1,6 +1,7 @@
 test_that("Simulating with table of parameters works", {
   skip_on_cran()
   parameters_table <- data.frame(
+    KA = rnorm(10, 1, .1),
     CL = rnorm(10, 5, 5),
     V = rnorm(10, 5, 0.5)
   )

From 3b024308dc09e7aa0cf380d6e79ce9e908af7406 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 11:31:11 -0800
Subject: [PATCH 07/10] move more tests

---
 R/sim.R                                       | 100 ---
 tests/testthat/test_bioav_def.R               |  44 --
 ...rameters.R => test_calculate_parameters.R} |   0
 tests/testthat/test_multi_obs.R               | 229 -------
 tests/testthat/test_multiple_regimens.R       |  24 -
 tests/testthat/test_new_ode_model.R           |  43 ++
 ...translation.R => test_pkpdsim_to_nlmixr.R} |   0
 tests/testthat/test_sim.R                     | 568 +++++++++++++++++-
 tests/testthat/test_sim_lagtime.R             |  26 -
 tests/testthat/test_t_init.R                  |  31 -
 tests/testthat/test_time_rounding_bug.R       |  15 -
 tests/testthat/test_timevar_cov.R             | 111 ----
 12 files changed, 610 insertions(+), 581 deletions(-)
 delete mode 100644 tests/testthat/test_bioav_def.R
 rename tests/testthat/{test_calc_parameters.R => test_calculate_parameters.R} (100%)
 delete mode 100644 tests/testthat/test_multi_obs.R
 delete mode 100644 tests/testthat/test_multiple_regimens.R
 rename tests/testthat/{test_nlmixr_translation.R => test_pkpdsim_to_nlmixr.R} (100%)
 delete mode 100644 tests/testthat/test_sim_lagtime.R
 delete mode 100644 tests/testthat/test_t_init.R
 delete mode 100644 tests/testthat/test_time_rounding_bug.R
 delete mode 100644 tests/testthat/test_timevar_cov.R

diff --git a/R/sim.R b/R/sim.R
index b8eb9a26..feda5c67 100644
--- a/R/sim.R
+++ b/R/sim.R
@@ -653,103 +653,3 @@ sim <- function (ode = NULL,
   attr(comb, "parameters") <- attr(ode, "parameters")
   return(comb)
 }
-
-describe("Mixture models", {
-  # Uses model and covariates defined in setup.R:
-  # - mod_mixture
-  # - covs_mixture
-
-  par_mixture <- list(CL = 3, V = 50)
-  reg_mixture <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
-  t_obs_mixture <- seq(0, 36, 4)
-
-  test_that("mixture model works properly for single patient", {
-    res0 <- sim_ode(mod_mixture, parameters = par_mixture, regimen = reg_mixture, covariates = covs_mixture, t_obs = t_obs_mixture, only_obs=T) # mixture_group not supplied
-    res1 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T)
-    res2 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 2, only_obs=T)
-    expect_equal(round(res0[res0$t == 24,]$y, 2), 9.07) # should use whatever is in `parameters`
-    expect_equal(round(res1[res1$t == 24,]$y, 2), 5.82)
-    expect_equal(round(res2[res2$t == 24,]$y, 2), 1.15)
-  })
-
-  test_that("mixture model works properly when vectorized (using parameters_table)", {
-    partab <- data.frame(CL = rep(0, 6), V = rep(50, 6))
-    suppressMessages({
-      expect_error(sim_ode(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T))
-      res1 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(1, 6), only_obs=T)
-      res2 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
-    })
-    expect_equal(round(res1[res1$t == 24,]$y, 2), rep(5.82, 6))
-    expect_equal(round(res2[res2$t == 24,]$y, 2), rep(c(5.82, 1.15), 3))
-    expect_equal(res2[res2$id == 1,]$CL[1], 5)
-    expect_equal(res2[res2$id == 2,]$CL[1], 15)
-    expect_equal(res2[res2$id == 3,]$CL[1], 5)
-  })
-
-  test_that("mixture model works properly when vectorized (using covariates_table)", {
-    covtab <- data.frame(ID = 1:8, WT = rep(seq(40, 130, 30), 2))
-    suppressMessages({
-      expect_error(sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = 1, only_obs=T))
-      res <- sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = rep(c(1, 2), each=4), only_obs=T)
-    })
-    expect_equal(round(res[res$t == 24,]$y, 2), c(9.39, 5.82, 3.83, 2.65, 2.99, 1.15, 0.52, 0.25))
-  })
-
-})
-
-describe("Compartment mapping", {
-  # Uses model defined in setup.R:
-  # - pk1cmt_oral_cmt_mapping
-
-  test_that("Compartment mapping is added to attributes", {
-    expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["oral"]], 1)
-    expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["infusion"]], 2)
-  })
-
-  test_that("Admin route is interpreted and simulated correctly", {
-    regimen <- new_regimen(
-      amt = c(100, 100, 100, 100),
-      times = c(0, 12, 24, 36),
-      type = c("oral", "oral", "infusion", "infusion"),
-      t_inf = c(0, 0, 1, 1)
-    )
-    p <- list(KA = 1, CL = 5, V = 50)
-    res <-  sim_ode(
-      ode = pk1cmt_oral_cmt_mapping,
-      parameters = p,
-      regimen = regimen,
-      only_obs = FALSE
-    )
-    expect_equal(round(res$y[res$comp == 1 & res$t == 25], 4), 2e-04)
-    expect_true(res$y[res$comp == 2 & res$t == 25] >= 100)
-  })
-
-  test_that("multiple scaling types on one compartment works", {
-    skip_on_cran()
-    mod <- new_ode_model(
-      code = "
-        dAdt[1] = -KA * A[1];
-        dAdt[2] = -(CL/V) * A[2] + KA*A[1];
-      ",
-      obs = list(
-        cmt = c(2, 2),
-        scale = c(1, "V"),
-        label = c("abs", "conc")
-      ),
-      cpp_show_code = FALSE
-    )
-    par <- list(CL = 5, V = 50, KA = .5)
-    reg <- new_regimen(amt = 100, n = 5, interval = 12)
-    res <- sim_ode(
-      ode = mod,
-      parameters = par,
-      regimen = reg,
-      only_obs = TRUE
-    )
-    dat <- cbind(res[res$comp == "abs",]$y, res[res$comp == "conc",]$y)
-    expect_true("PKPDsim" %in% class(mod))
-    expect_equal(length(unique(res$comp)), 2)
-    expect_equal(round(dat[,1],1), round(dat[,2]*par$V,1))
-  })
-
-})
\ No newline at end of file
diff --git a/tests/testthat/test_bioav_def.R b/tests/testthat/test_bioav_def.R
deleted file mode 100644
index c4d7674a..00000000
--- a/tests/testthat/test_bioav_def.R
+++ /dev/null
@@ -1,44 +0,0 @@
-parameters <- list(KA = 0.5, CL = 5, V = 50)
-covs <- list(WT = new_covariate(50))
-reg <- new_regimen(amt = 100, n = 1, interval = 12, type = "bolus")
-mod1 <- oral_1cmt_allometric # defined in setup.R using same covs and pars as above
-y1 <- sim_ode(mod1, parameters = parameters, regimen = reg, covariates = covs, only_obs=TRUE)$y
-
-test_that("bioav numeric option working", {
-  skip_on_cran()
-  mod2 <- new_ode_model(
-    code = "
-    CLi = CL * pow(WT/70, 0.75)
-    dAdt[1] = -KA * A[1]
-    dAdt[2] = KA*A[1] - (CLi/V)*A[2]
-  ",
-  dose = list(cmt = 1, bioav = 0.5),
-  covariates = covs,
-  declare_variables = "CLi",
-  parameters = parameters
-  )
-
-  y2 <- sim_ode(mod2, parameters = parameters, regimen = reg, covariates = covs, only_obs=TRUE)$y
-
-  expect_equal(round(y1,1),  round(y2*2, 1))
-})
-
-test_that("bioav string option working", {
-  skip_on_cran()
-  mod3 <- new_ode_model(
-    code = "
-    CLi = CL * pow(WT/70, 0.75)
-    dAdt[1] = -KA * A[1]
-    dAdt[2] = KA*A[1] - (CLi/V)*A[2]
-  ",
-  dose = list(cmt = 1, bioav = "WT/70"),
-  covariates = covs,
-  declare_variables = "CLi",
-  parameters = parameters
-  )
-
-  y3 <- sim_ode(mod3, parameters = parameters, regimen = reg, covariates = covs, only_obs=TRUE)$y
-
-  expect_equal(round(y1*(50/70),1),  round(y3, 1))
-})
-
diff --git a/tests/testthat/test_calc_parameters.R b/tests/testthat/test_calculate_parameters.R
similarity index 100%
rename from tests/testthat/test_calc_parameters.R
rename to tests/testthat/test_calculate_parameters.R
diff --git a/tests/testthat/test_multi_obs.R b/tests/testthat/test_multi_obs.R
deleted file mode 100644
index 15be1909..00000000
--- a/tests/testthat/test_multi_obs.R
+++ /dev/null
@@ -1,229 +0,0 @@
-# Uses model defined in setup.R:
-# - pk_multi_obs
-# - vars_multi_obs
-
-regimen_multi_obs <- new_regimen(amt = 100, interval = 12, n = 5, type="infusion", t_inf = 1)
-parameters_multi_obs <- list("CL" = 15, "V" = 150)
-omega_multi_obs <- PKPDsim::cv_to_omega(list("CL" = 0.2, "V" = 0.2), parameters_multi_obs[1:2])
-
-test_that("obs types are output by `sim`", {
-  obs_type <- c(1,2,1,3,1)
-  data <- sim(
-    ode = pk_multi_obs,
-    parameters = list(CL = 20, V = 200),
-    regimen = regimen_multi_obs,
-    int_step_size = 0.1,
-    only_obs = TRUE,
-    obs_type = obs_type,
-    t_obs = c(2, 4, 6, 8, 12),
-    output_include = list("variables" = TRUE)
-  )
-  expect_equal(data$obs_type, obs_type)
-  expect_equal(data$y, diag(as.matrix(data[1:5,5+obs_type])))
-})
-
-
-test_that("check obs at same timepoint but with different obs_type", {
-  obs_type <- c(1,2,1,3,1)
-  t_same <- sim(
-    ode = pk_multi_obs,
-    parameters = list(CL = 20, V = 200),
-    regimen = regimen_multi_obs,
-    int_step_size = 0.1,
-    only_obs = TRUE,
-    obs_type = obs_type,
-    t_obs = c(2, 4, 4, 8, 8),
-    output_include = list("variables" = TRUE)
-  )
-  expect_equal(t_same$t[2], t_same$t[3])
-  expect_equal(t_same$t[4], t_same$t[5])
-  expect_equal(t_same$y[3], t_same$METAB[3])
-  expect_equal(t_same$y[2], t_same$CONC[2])
-  expect_equal(t_same$y[5], t_same$METAB2[5])
-  expect_equal(t_same$y[4], t_same$CONC[4])
-})
-
-
-test_that("check that residual error correctly applied to right var", {
-  obs_type <- c(1,2,1,3,1)
-  set.seed(12345)
-  ruv_term3 <- list(prop = c(0, 0, 0.1), add = c(0, 0, 0.1))
-  ruv_term1 <- list(prop = c(.1, 0, 0), add = c(1, 0, 0))
-
-  data2 <- sim(
-    ode = pk_multi_obs,
-    parameters = list(CL = 20, V = 200),
-    regimen = regimen_multi_obs,
-    int_step_size = 0.1,
-    only_obs = TRUE,
-    obs_type = obs_type,
-    t_obs = c(2, 4, 6, 8, 12),
-    output_include = list("variables" = TRUE),
-    res_var = ruv_term3
-  )
-  y <- diag(as.matrix(data2[1:5, 5 + obs_type]))
-
-  # no residual error for obs_type 1 and 2, only 3
-  expect_equal(data2$y[-4], y[-4])
-  expect_true(data2$y[-4][4] != y[4])
-
-  data3 <- sim(
-    ode = pk_multi_obs,
-    parameters = list(CL = 20, V = 200),
-    regimen = regimen_multi_obs,
-    int_step_size = 0.1,
-    only_obs = TRUE,
-    obs_type = obs_type,
-    t_obs = c(2, 4, 6, 8, 12),
-    output_include = list("variables" = TRUE),
-    res_var = ruv_term1
-  )
-  y <- diag(as.matrix(data2[1:5, 5 + obs_type]))
-
-  # no residual error for obs_type 2/3, only and 1
-  expect_equal(data3$y[-c(1,3,5)], y[-c(1,3,5)])
-  expect_true(all(data3$y[c(1,3,5)] != y[c(1,3,5)]))
-})
-
-
-test_that("specifying ruv as multi-type when only 1 obs_type", {
-  obs_type <- c(1,2,1,3,1)
-  ruv_single <- list(prop = 0.1, add = 1)
-  ruv_multi <- list(prop = c(0.1, 1), add = c(1, 20))
-
-  s_single <- sim(
-    ode = pk_multi_obs,
-    parameters = parameters_multi_obs,
-    n_ind = 1,
-    regimen = regimen_multi_obs,
-    only_obs = TRUE,
-    t_obs = c(2,4,6,8),
-    seed = 123456,
-    res_var = ruv_single
-  )
-
-  ## specified as multi, but obs_type is all 1, so should
-  ## give same results as s_single
-  s_multi1 <- sim(
-    ode = pk_multi_obs,
-    parameters = parameters_multi_obs,
-    n_ind = 1,
-    regimen = regimen_multi_obs,
-    only_obs = TRUE,
-    t_obs = c(2,4,6,8),
-    obs_type = c(1, 1, 1, 1),
-    seed = 123456,
-    res_var = ruv_multi,
-    t_init = 10
-  )
-  expect_equal(s_multi1$y, s_single$y)
-
-  ## specifying ruv as multi-type when multiple obs_type
-  ## should now give different results
-  s_multi1 <- sim(
-    ode = pk_multi_obs,
-    parameters = parameters_multi_obs,
-    n_ind = 1,
-    regimen = regimen_multi_obs,
-    only_obs = TRUE,
-    t_obs = c(2,4,6,8),
-    obs_type = c(1, 2, 1, 2),
-    seed = 123456,
-    res_var = ruv_multi,
-    t_init = 10
-  )
-  expect_equal(s_multi1$y[c(1, 3)], s_single$y[c(1,3)])
-  expect_true(sum(abs(s_single$y[c(2,4)] - s_multi1$y[c(2,4)])) > 50)
-})
-
-test_that("multi-obs with baseline and obs_time = dose_time works correctly", {
-  tmp <- sim(
-    pk_multi_obs,
-    parameters = parameters_multi_obs,
-    regimen = regimen_multi_obs,
-    t_obs = c(0, 0, 6, 6),
-    obs_type = c(1, 4, 1, 4),
-    only_obs = TRUE,
-    output_include = list(variables = TRUE),
-    return_design = F
-  )
-  expect_equal(
-    tmp$y[tmp$obs_type == 1],
-    tmp$CONC[tmp$obs_type == 1]
-  )
-  expect_equal(
-    tmp$y[tmp$obs_type == 4],
-    tmp$ACT[tmp$obs_type == 4]
-  )
-})
-
-test_that("multi-obs model simulates correctly (bug 202205)", {
-  skip_on_cran()
-
-  pars <- list(
-    CL = 23.9,
-    CLM = 5.19,
-    V = 107,
-    Q = 3.31,
-    V2 = 46.2,
-    VM = 111,
-    KA = 18.6,
-    TLAG = 0.415
-  )
-  covs <- list(WT = PKPDsim::new_covariate(80))
-  mod <- new_ode_model(
-    code = "dAdt[0] = -KAi*A[0] \
-    dAdt[1] = +KAi*A[0] - (CLi/Vi)*A[1] - (Qi/Vi)*A[1] + (Qi/V2i)*A[2] \
-    dAdt[2] = +(Qi/Vi)*A[1] - (Qi/V2i)*A[2] \
-    dAdt[3] = +(CLi/Vi)*A[1] - (CLMi/VMi)*A[3] \
-    dAdt[4] = A[1]*1000.0/Vi \
-    CONC = A[1]/(Vi/1000.0) \
-    CONCM = A[3]/(VMi/1000.0) \
-    CONCT = CONC+CONCM \
-  ",
-    pk = "KAi = KA \
-    CLi = CL * pow(WT/70.0, 0.75) \
-    Vi = V *(WT/70.0) \
-    Qi = Q * pow(WT/70.0, 0.75)\
-    V2i = V2 * (WT/70.0) \
-    CLMi = CLM * pow(WT/70.0, 0.75) \
-    VMi = VM *(WT/70.0) \
-  ",
-    parameters = pars,
-    declare_variables = c( "CLi", "Vi", "V2i", "Qi", "KAi", "VMi", "CLMi", "CONC", "CONCM", "CONCT" ),
-    obs = list(
-      variable = c( "CONC", "CONCM", "CONC", "CONCM" )
-    ),
-    covariates = covs
-  )
-  regimen <- new_regimen(amt = 3.7, n=10, interval = 24, type = "oral", cmt = 1)
-  data <- data.frame(
-    t = 70,
-    y = 1,
-    evid = 0,
-    loq = 0,
-    obs_type = 2,
-    dv = 1
-  )
-  t_obs <- seq(0, 50, .5)
-  pop <- sim(
-    ode = mod,
-    regimen = regimen,
-    parameters = pars,
-    covariates = covs,
-    t_obs = rep(t_obs, each = 4),
-    only_obs = TRUE,
-    obs_type = rep(1:4, length(t_obs))
-  )
-  # correct obs_types
-  expect_equal(
-    pop[pop$t %in% c(23, 47), ]$obs_type,
-    c(1,2,3,4,1,2,3,4)
-  )
-  # correct output values
-  expect_equal(
-    round(pop[pop$t %in% c(23, 47), ]$y, 2),
-    c(0.52, 12.34, 0.52, 12.34, 0.63, 17.17, 0.63, 17.17)
-  )
-
-})
diff --git a/tests/testthat/test_multiple_regimens.R b/tests/testthat/test_multiple_regimens.R
deleted file mode 100644
index 59a312b1..00000000
--- a/tests/testthat/test_multiple_regimens.R
+++ /dev/null
@@ -1,24 +0,0 @@
-test_that("multiple regimens for multiple individuals can be simulated", {
-  # Uses a model defined in `setup.R`
-  cov_table <- data.frame(WT = rnorm(10, 70, 5))
-  multi_regs <- list()
-  for(i in seq(cov_table$WT)) {
-    multi_regs[[i]] <- new_regimen(amt = 10 * cov_table$WT[i], interval = 12, type = "infusion")
-  }
-  class(multi_regs) <- "regimen_multiple"
-
-  par <- list(CL = 5, V = 50)
-  reg <- new_regimen(amt = 2000, interval = 24, type = "infusion")
-  covariates = list(WT = new_covariate(1))
-  res <- sim_ode(
-    ode = mod_1cmt_iv,
-    parameters = par,
-    covariates = covariates,
-    regimen = multi_regs,
-    only_obs = TRUE
-  )
-  expect_equal(length(unique(res$id)),10)
-  expect_equal(sum(is.na(res$y)), 0)
-  # IDs ordered correctly:
-  expect_true(length(unique(res$id)) == 10 && all(diff(res$id) >= 0))
-})
diff --git a/tests/testthat/test_new_ode_model.R b/tests/testthat/test_new_ode_model.R
index c27f3205..0b3c101f 100644
--- a/tests/testthat/test_new_ode_model.R
+++ b/tests/testthat/test_new_ode_model.R
@@ -225,3 +225,46 @@ test_that("IOV: Change in F in 2nd bin is applied in 2nd bin and not later.", {
   expect_true(min(res1[res1$y != res2$y,]$t) <= 25)
 })
 
+describe("Models with bioavailability", {
+  parameters <- list(KA = 0.5, CL = 5, V = 50)
+  covs <- list(WT = new_covariate(50))
+  reg <- new_regimen(amt = 100, n = 1, interval = 12, type = "bolus")
+  mod1 <- oral_1cmt_allometric # defined in setup.R using same covs and pars as above
+  y1 <- sim_ode(mod1, parameters = parameters, regimen = reg, covariates = covs, only_obs=TRUE)$y
+
+  test_that("bioav numeric option working", {
+    skip_on_cran()
+    mod2 <- new_ode_model(
+      code = "
+      CLi = CL * pow(WT/70, 0.75)
+      dAdt[1] = -KA * A[1]
+      dAdt[2] = KA*A[1] - (CLi/V)*A[2]
+    ",
+    dose = list(cmt = 1, bioav = 0.5),
+    covariates = covs,
+    declare_variables = "CLi",
+    parameters = parameters
+    )
+
+    y2 <- sim_ode(mod2, parameters = parameters, regimen = reg, covariates = covs, only_obs=TRUE)$y
+
+    expect_equal(round(y1,1),  round(y2*2, 1))
+  })
+
+  test_that("bioav string option working", {
+    skip_on_cran()
+    mod3 <- new_ode_model(
+      code = "
+      CLi = CL * pow(WT/70, 0.75)
+      dAdt[1] = -KA * A[1]
+      dAdt[2] = KA*A[1] - (CLi/V)*A[2]
+    ",
+    dose = list(cmt = 1, bioav = "WT/70"),
+    covariates = covs,
+    declare_variables = "CLi",
+    parameters = parameters
+    )
+    y3 <- sim_ode(mod3, parameters = parameters, regimen = reg, covariates = covs, only_obs=TRUE)$y
+    expect_equal(round(y1*(50/70),1),  round(y3, 1))
+  })
+})
\ No newline at end of file
diff --git a/tests/testthat/test_nlmixr_translation.R b/tests/testthat/test_pkpdsim_to_nlmixr.R
similarity index 100%
rename from tests/testthat/test_nlmixr_translation.R
rename to tests/testthat/test_pkpdsim_to_nlmixr.R
diff --git a/tests/testthat/test_sim.R b/tests/testthat/test_sim.R
index 9abb3789..d8984bec 100644
--- a/tests/testthat/test_sim.R
+++ b/tests/testthat/test_sim.R
@@ -743,7 +743,7 @@ describe("IOV", {
 })
 
 describe("Compare results from sims with references", {
-  
+
   # Uses models defined in setup.R:
   # - mod_1cmt_oral
   # - mod_1cmt_iv
@@ -999,4 +999,570 @@ describe("Compare results from sims with references", {
     expect_equal(length(unique(dat$id)), 3)
   })
 
+})
+
+describe("Simulate with multiple observation types", {
+
+  # Uses model defined in setup.R:
+  # - pk_multi_obs
+  # - vars_multi_obs
+
+  regimen_multi_obs <- new_regimen(amt = 100, interval = 12, n = 5, type="infusion", t_inf = 1)
+  parameters_multi_obs <- list("CL" = 15, "V" = 150)
+  omega_multi_obs <- PKPDsim::cv_to_omega(list("CL" = 0.2, "V" = 0.2), parameters_multi_obs[1:2])
+
+  test_that("obs types are output by `sim`", {
+    obs_type <- c(1,2,1,3,1)
+    data <- sim(
+      ode = pk_multi_obs,
+      parameters = list(CL = 20, V = 200),
+      regimen = regimen_multi_obs,
+      int_step_size = 0.1,
+      only_obs = TRUE,
+      obs_type = obs_type,
+      t_obs = c(2, 4, 6, 8, 12),
+      output_include = list("variables" = TRUE)
+    )
+    expect_equal(data$obs_type, obs_type)
+    expect_equal(data$y, diag(as.matrix(data[1:5,5+obs_type])))
+  })
+
+
+  test_that("check obs at same timepoint but with different obs_type", {
+    obs_type <- c(1,2,1,3,1)
+    t_same <- sim(
+      ode = pk_multi_obs,
+      parameters = list(CL = 20, V = 200),
+      regimen = regimen_multi_obs,
+      int_step_size = 0.1,
+      only_obs = TRUE,
+      obs_type = obs_type,
+      t_obs = c(2, 4, 4, 8, 8),
+      output_include = list("variables" = TRUE)
+    )
+    expect_equal(t_same$t[2], t_same$t[3])
+    expect_equal(t_same$t[4], t_same$t[5])
+    expect_equal(t_same$y[3], t_same$METAB[3])
+    expect_equal(t_same$y[2], t_same$CONC[2])
+    expect_equal(t_same$y[5], t_same$METAB2[5])
+    expect_equal(t_same$y[4], t_same$CONC[4])
+  })
+
+
+  test_that("check that residual error correctly applied to right var", {
+    obs_type <- c(1,2,1,3,1)
+    set.seed(12345)
+    ruv_term3 <- list(prop = c(0, 0, 0.1), add = c(0, 0, 0.1))
+    ruv_term1 <- list(prop = c(.1, 0, 0), add = c(1, 0, 0))
+
+    data2 <- sim(
+      ode = pk_multi_obs,
+      parameters = list(CL = 20, V = 200),
+      regimen = regimen_multi_obs,
+      int_step_size = 0.1,
+      only_obs = TRUE,
+      obs_type = obs_type,
+      t_obs = c(2, 4, 6, 8, 12),
+      output_include = list("variables" = TRUE),
+      res_var = ruv_term3
+    )
+    y <- diag(as.matrix(data2[1:5, 5 + obs_type]))
+
+    # no residual error for obs_type 1 and 2, only 3
+    expect_equal(data2$y[-4], y[-4])
+    expect_true(data2$y[-4][4] != y[4])
+
+    data3 <- sim(
+      ode = pk_multi_obs,
+      parameters = list(CL = 20, V = 200),
+      regimen = regimen_multi_obs,
+      int_step_size = 0.1,
+      only_obs = TRUE,
+      obs_type = obs_type,
+      t_obs = c(2, 4, 6, 8, 12),
+      output_include = list("variables" = TRUE),
+      res_var = ruv_term1
+    )
+    y <- diag(as.matrix(data2[1:5, 5 + obs_type]))
+
+    # no residual error for obs_type 2/3, only and 1
+    expect_equal(data3$y[-c(1,3,5)], y[-c(1,3,5)])
+    expect_true(all(data3$y[c(1,3,5)] != y[c(1,3,5)]))
+  })
+
+
+  test_that("specifying ruv as multi-type when only 1 obs_type", {
+    obs_type <- c(1,2,1,3,1)
+    ruv_single <- list(prop = 0.1, add = 1)
+    ruv_multi <- list(prop = c(0.1, 1), add = c(1, 20))
+
+    s_single <- sim(
+      ode = pk_multi_obs,
+      parameters = parameters_multi_obs,
+      n_ind = 1,
+      regimen = regimen_multi_obs,
+      only_obs = TRUE,
+      t_obs = c(2,4,6,8),
+      seed = 123456,
+      res_var = ruv_single
+    )
+
+    ## specified as multi, but obs_type is all 1, so should
+    ## give same results as s_single
+    s_multi1 <- sim(
+      ode = pk_multi_obs,
+      parameters = parameters_multi_obs,
+      n_ind = 1,
+      regimen = regimen_multi_obs,
+      only_obs = TRUE,
+      t_obs = c(2,4,6,8),
+      obs_type = c(1, 1, 1, 1),
+      seed = 123456,
+      res_var = ruv_multi,
+      t_init = 10
+    )
+    expect_equal(s_multi1$y, s_single$y)
+
+    ## specifying ruv as multi-type when multiple obs_type
+    ## should now give different results
+    s_multi1 <- sim(
+      ode = pk_multi_obs,
+      parameters = parameters_multi_obs,
+      n_ind = 1,
+      regimen = regimen_multi_obs,
+      only_obs = TRUE,
+      t_obs = c(2,4,6,8),
+      obs_type = c(1, 2, 1, 2),
+      seed = 123456,
+      res_var = ruv_multi,
+      t_init = 10
+    )
+    expect_equal(s_multi1$y[c(1, 3)], s_single$y[c(1,3)])
+    expect_true(sum(abs(s_single$y[c(2,4)] - s_multi1$y[c(2,4)])) > 50)
+  })
+
+  test_that("multi-obs with baseline and obs_time = dose_time works correctly", {
+    tmp <- sim(
+      pk_multi_obs,
+      parameters = parameters_multi_obs,
+      regimen = regimen_multi_obs,
+      t_obs = c(0, 0, 6, 6),
+      obs_type = c(1, 4, 1, 4),
+      only_obs = TRUE,
+      output_include = list(variables = TRUE),
+      return_design = F
+    )
+    expect_equal(
+      tmp$y[tmp$obs_type == 1],
+      tmp$CONC[tmp$obs_type == 1]
+    )
+    expect_equal(
+      tmp$y[tmp$obs_type == 4],
+      tmp$ACT[tmp$obs_type == 4]
+    )
+  })
+
+  test_that("multi-obs model simulates correctly (bug 202205)", {
+    skip_on_cran()
+
+    pars <- list(
+      CL = 23.9,
+      CLM = 5.19,
+      V = 107,
+      Q = 3.31,
+      V2 = 46.2,
+      VM = 111,
+      KA = 18.6,
+      TLAG = 0.415
+    )
+    covs <- list(WT = PKPDsim::new_covariate(80))
+    mod <- new_ode_model(
+      code = "dAdt[0] = -KAi*A[0] \
+      dAdt[1] = +KAi*A[0] - (CLi/Vi)*A[1] - (Qi/Vi)*A[1] + (Qi/V2i)*A[2] \
+      dAdt[2] = +(Qi/Vi)*A[1] - (Qi/V2i)*A[2] \
+      dAdt[3] = +(CLi/Vi)*A[1] - (CLMi/VMi)*A[3] \
+      dAdt[4] = A[1]*1000.0/Vi \
+      CONC = A[1]/(Vi/1000.0) \
+      CONCM = A[3]/(VMi/1000.0) \
+      CONCT = CONC+CONCM \
+    ",
+      pk = "KAi = KA \
+      CLi = CL * pow(WT/70.0, 0.75) \
+      Vi = V *(WT/70.0) \
+      Qi = Q * pow(WT/70.0, 0.75)\
+      V2i = V2 * (WT/70.0) \
+      CLMi = CLM * pow(WT/70.0, 0.75) \
+      VMi = VM *(WT/70.0) \
+    ",
+      parameters = pars,
+      declare_variables = c( "CLi", "Vi", "V2i", "Qi", "KAi", "VMi", "CLMi", "CONC", "CONCM", "CONCT" ),
+      obs = list(
+        variable = c( "CONC", "CONCM", "CONC", "CONCM" )
+      ),
+      covariates = covs
+    )
+    regimen <- new_regimen(amt = 3.7, n=10, interval = 24, type = "oral", cmt = 1)
+    data <- data.frame(
+      t = 70,
+      y = 1,
+      evid = 0,
+      loq = 0,
+      obs_type = 2,
+      dv = 1
+    )
+    t_obs <- seq(0, 50, .5)
+    pop <- sim(
+      ode = mod,
+      regimen = regimen,
+      parameters = pars,
+      covariates = covs,
+      t_obs = rep(t_obs, each = 4),
+      only_obs = TRUE,
+      obs_type = rep(1:4, length(t_obs))
+    )
+    # correct obs_types
+    expect_equal(
+      pop[pop$t %in% c(23, 47), ]$obs_type,
+      c(1,2,3,4,1,2,3,4)
+    )
+    # correct output values
+    expect_equal(
+      round(pop[pop$t %in% c(23, 47), ]$y, 2),
+      c(0.52, 12.34, 0.52, 12.34, 0.63, 17.17, 0.63, 17.17)
+    )
+
+  })
+
+})
+
+describe("Multiple regimens", {
+  test_that("multiple regimens for multiple individuals can be simulated", {
+
+    # Uses a model defined in `setup.R`
+    cov_table <- data.frame(WT = rnorm(10, 70, 5))
+    multi_regs <- list()
+    for(i in seq(cov_table$WT)) {
+      multi_regs[[i]] <- new_regimen(amt = 10 * cov_table$WT[i], interval = 12, type = "infusion")
+    }
+    class(multi_regs) <- "regimen_multiple"
+
+    par <- list(CL = 5, V = 50)
+    reg <- new_regimen(amt = 2000, interval = 24, type = "infusion")
+    covariates = list(WT = new_covariate(1))
+    res <- sim_ode(
+      ode = mod_1cmt_iv,
+      parameters = par,
+      covariates = covariates,
+      regimen = multi_regs,
+      only_obs = TRUE
+    )
+    expect_equal(length(unique(res$id)),10)
+    expect_equal(sum(is.na(res$y)), 0)
+    # IDs ordered correctly:
+    expect_true(length(unique(res$id)) == 10 && all(diff(res$id) >= 0))
+  })
+
+})
+
+
+describe("Mixture models", {
+  # Uses model and covariates defined in setup.R:
+  # - mod_mixture
+  # - covs_mixture
+
+  par_mixture <- list(CL = 3, V = 50)
+  reg_mixture <- new_regimen(amt = 250, n = 5, interval = 6, type = 'infusion', t_inf = 1)
+  t_obs_mixture <- seq(0, 36, 4)
+
+  test_that("mixture model works properly for single patient", {
+    res0 <- sim_ode(mod_mixture, parameters = par_mixture, regimen = reg_mixture, covariates = covs_mixture, t_obs = t_obs_mixture, only_obs=T) # mixture_group not supplied
+    res1 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T)
+    res2 <- sim(mod_mixture, parameters = par_mixture, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 2, only_obs=T)
+    expect_equal(round(res0[res0$t == 24,]$y, 2), 9.07) # should use whatever is in `parameters`
+    expect_equal(round(res1[res1$t == 24,]$y, 2), 5.82)
+    expect_equal(round(res2[res2$t == 24,]$y, 2), 1.15)
+  })
+
+  test_that("mixture model works properly when vectorized (using parameters_table)", {
+    partab <- data.frame(CL = rep(0, 6), V = rep(50, 6))
+    suppressMessages({
+      expect_error(sim_ode(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = 1, only_obs=T))
+      res1 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(1, 6), only_obs=T)
+      res2 <- sim(mod_mixture, parameters_table = partab, regimen = reg_mixture, t_obs = t_obs_mixture, covariates = covs_mixture, mixture_group = rep(c(1,2), 3), only_obs=T, output_include = list(parameters = TRUE))
+    })
+    expect_equal(round(res1[res1$t == 24,]$y, 2), rep(5.82, 6))
+    expect_equal(round(res2[res2$t == 24,]$y, 2), rep(c(5.82, 1.15), 3))
+    expect_equal(res2[res2$id == 1,]$CL[1], 5)
+    expect_equal(res2[res2$id == 2,]$CL[1], 15)
+    expect_equal(res2[res2$id == 3,]$CL[1], 5)
+  })
+
+  test_that("mixture model works properly when vectorized (using covariates_table)", {
+    covtab <- data.frame(ID = 1:8, WT = rep(seq(40, 130, 30), 2))
+    suppressMessages({
+      expect_error(sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = 1, only_obs=T))
+      res <- sim(mod_mixture, parameters = par_mixture, covariates_table = covtab, regimen = reg_mixture, t_obs = t_obs_mixture, mixture_group = rep(c(1, 2), each=4), only_obs=T)
+    })
+    expect_equal(round(res[res$t == 24,]$y, 2), c(9.39, 5.82, 3.83, 2.65, 2.99, 1.15, 0.52, 0.25))
+  })
+
+})
+
+describe("Compartment mapping", {
+  # Uses model defined in setup.R:
+  # - pk1cmt_oral_cmt_mapping
+
+  test_that("Compartment mapping is added to attributes", {
+    expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["oral"]], 1)
+    expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["infusion"]], 2)
+  })
+
+  test_that("Admin route is interpreted and simulated correctly", {
+    regimen <- new_regimen(
+      amt = c(100, 100, 100, 100),
+      times = c(0, 12, 24, 36),
+      type = c("oral", "oral", "infusion", "infusion"),
+      t_inf = c(0, 0, 1, 1)
+    )
+    p <- list(KA = 1, CL = 5, V = 50)
+    res <-  sim_ode(
+      ode = pk1cmt_oral_cmt_mapping,
+      parameters = p,
+      regimen = regimen,
+      only_obs = FALSE
+    )
+    expect_equal(round(res$y[res$comp == 1 & res$t == 25], 4), 2e-04)
+    expect_true(res$y[res$comp == 2 & res$t == 25] >= 100)
+  })
+
+  test_that("multiple scaling types on one compartment works", {
+    skip_on_cran()
+    mod <- new_ode_model(
+      code = "
+        dAdt[1] = -KA * A[1];
+        dAdt[2] = -(CL/V) * A[2] + KA*A[1];
+      ",
+      obs = list(
+        cmt = c(2, 2),
+        scale = c(1, "V"),
+        label = c("abs", "conc")
+      ),
+      cpp_show_code = FALSE
+    )
+    par <- list(CL = 5, V = 50, KA = .5)
+    reg <- new_regimen(amt = 100, n = 5, interval = 12)
+    res <- sim_ode(
+      ode = mod,
+      parameters = par,
+      regimen = reg,
+      only_obs = TRUE
+    )
+    dat <- cbind(res[res$comp == "abs",]$y, res[res$comp == "conc",]$y)
+    expect_true("PKPDsim" %in% class(mod))
+    expect_equal(length(unique(res$comp)), 2)
+    expect_equal(round(dat[,1],1), round(dat[,2]*par$V,1))
+  })
+
+})
+
+describe("Model simulation with lagtime", {
+
+  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))
+  })
+
+})
+
+
+describe("Models with t_init", {
+  # Uses model defined in setup.R (conditional, NOT_CRAN only):
+  # - mod_t_init
+
+  skip_on_cran()
+
+  # Test t_init functionality
+  ## e.g. TDM before first dose:
+  ## at t=-8, conc=10000
+  ## Use this as true value in the simulations
+  par_t_init <- list(CL = 7.67, V = 97.7, TDM_INIT = 500)
+  reg_t_init <- new_regimen(amt = 100000, times=c(0, 24), type="bolus")
+  s <- sim(
+    ode = mod_t_init,
+    parameters = par_t_init,
+    n_ind = 1,
+    regimen = reg_t_init,
+    only_obs = TRUE,
+    output_include = list(variables = TRUE),
+    t_init = 10
+  )
+
+  test_that("TDM before first dose is considered a true initial value", {
+    expect_equal(sum(is.na(s$y)), 0)
+    expect_equal(s$y[1], 500)
+    expect_equal(round(s$y[3]), 427)
+    expect_equal(round(s$y[13]),1157)
+  })
+
+  test_that("Variables are set (also in first row) when TDM before first dose", {
+    expect_equal(round(s$CONC[1:5], 1), c(500, 462.2, 427.3, 395.1, 365.3))
+  })
+
+})
+
+describe("Time rounding issues", {
+
+  # rounding time should not produce NAs in sim
+  ## time-rounding bug 20170804
+
+  test_that("No NAs related to rounding", {
+    # Uses models defined in setup.R
+    p <- list(CL = 5, V = 50, Q = 10, V2 = 150)
+    r1 <- new_regimen(amt = 100, times = c(0, 24, 36), type = "infusion")
+    dat1 <- sim_ode(
+      ode = mod_2cmt_iv,
+      parameters = p,
+      regimen = r1,
+      t_obs=seq(0, 150, length.out = 100)
+    )
+    expect_equal(sum(is.na(dat1$y)), 0)
+  })
+
+})
+
+describe("Timevarying covariates are handled properly", {
+  # Uses model defined in setup.R (conditional, NOT_CRAN only):
+  # - mod_2cmt_timevar
+
+  skip_on_cran()
+
+  par_timevar <- list(CL = 3, V = 50, Q = 2.5, V2 = 70)
+  reg_timevar <- new_regimen(
+    amt = 250,
+    n = 60,
+    interval = 6,
+    type = 'infusion',
+    t_inf = 1
+  )
+
+  test_that("timevarying covariates handled", {
+    # CLi changes by several orders of magnitude after
+    # steady state is achieved, which should produce
+    # a new steady state that is much lower
+    covs <- list(
+      CRCL = new_covariate(
+        value = c(4.662744, 5.798767, 6.195943, 6.2, 600),
+        times = c(0, 18, 29, 209, 210),
+        implementation = "locf"
+      )
+    )
+    t_obs <- seq(0, 360, 0.1)
+    sim1 <- sim_ode(
+      mod_2cmt_timevar,
+      parameters = par_timevar,
+      regimen = reg_timevar,
+      covariates = covs,
+      only_obs = TRUE,
+      t_obs = t_obs,
+      output_include = list(parameters = TRUE, covariates = TRUE)
+    )
+    expect_equal(sim1$CRCL[sim1$t == 209], 6.2)
+    expect_equal(sim1$CRCL[sim1$t == 210], 600)
+    expect_true(sim1$y[sim1$t == 35 * 6] > 10 * sim1$y[sim1$t == 60 * 6])
+  })
+
+  test_that("timevarying covariates are interpolated and affect PK", {
+    covs_locf <- list(
+      CRCL = new_covariate(
+        times = c(0, 48, 96),
+        value = c(3, 10, 3),
+        implementation = "locf"
+      )
+    )
+    covs_inter <- list(
+      CRCL = new_covariate(
+        times = c(0, 48, 96),
+        value = c(3, 10, 3),
+        implementation = "interpolate"
+      )
+    )
+    t_obs <- seq(0, 120, 2)
+    sim2_inter <- sim_ode(
+      mod_2cmt_timevar,
+      parameters = par_timevar,
+      regimen = reg_timevar,
+      covariates = covs_inter,
+      only_obs = TRUE,
+      t_obs = t_obs,
+      output_include = list(parameters = TRUE, covariates = TRUE, variables = TRUE)
+    )
+    sim2_locf <- sim_ode(
+      mod_2cmt_timevar,
+      parameters = par_timevar,
+      regimen = reg_timevar,
+      covariates = covs_locf,
+      only_obs = TRUE,
+      t_obs = t_obs,
+      output_include = list(parameters = TRUE, covariates = TRUE, variables = TRUE)
+    )
+
+    ## Check covariate changing for inter, but not for locf
+    expect_equal(
+      round(sim2_inter$CRCL, 3)[1:10],
+      c(3, 3.292, 3.583, 3.875, 4.167, 4.458, 4.75, 5.042, 5.333, 5.625)
+    )
+    expect_equal(
+      round(sim2_locf$CRCL, 3)[1:10],
+      rep(3, 10)
+    )
+
+    ## Check interpolated covariates actually affect PK parameters
+    expect_equal(
+      round(sim2_inter$CLi, 3)[1:10],
+      c(6, 6.292, 6.583, 6.875, 7.167, 7.458, 7.75, 8.042, 8.333, 8.625)
+    )
+    expect_equal(
+      round(sim2_locf$CLi, 3)[1:10],
+      rep(6, 10)
+    )
+
+    ## Check interpolated covariates actually affect simulated conc
+    expect_equal(
+      round(sim2_inter$y, 3)[1:10],
+      c(0, 0.32, 0.611, 0.788, 1.205, 1.531, 1.708, 2.098, 2.379, 2.503)
+    )
+    expect_equal(
+      round(sim2_locf$y, 3)[1:10],
+      c(0, 0.321, 0.617, 0.805, 1.239, 1.598, 1.813, 2.251, 2.598, 2.792)
+    )
+
+    ## Visual check:
+    # ggplot(sim2_inter, aes(x = t, y = y)) +
+    #   geom_line() +
+    #   geom_line(data = sim2_locf, colour = "blue")
+
+  })
+
 })
\ No newline at end of file
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))
-})
diff --git a/tests/testthat/test_t_init.R b/tests/testthat/test_t_init.R
deleted file mode 100644
index abd4e2e5..00000000
--- a/tests/testthat/test_t_init.R
+++ /dev/null
@@ -1,31 +0,0 @@
-# Uses model defined in setup.R (conditional, NOT_CRAN only):
-# - mod_t_init
-
-skip_on_cran()
-
-# Test t_init functionality
-## e.g. TDM before first dose:
-## at t=-8, conc=10000
-## Use this as true value in the simulations
-par_t_init <- list(CL = 7.67, V = 97.7, TDM_INIT = 500)
-reg_t_init <- new_regimen(amt = 100000, times=c(0, 24), type="bolus")
-s <- sim(
-  ode = mod_t_init,
-  parameters = par_t_init,
-  n_ind = 1,
-  regimen = reg_t_init,
-  only_obs = TRUE,
-  output_include = list(variables = TRUE),
-  t_init = 10
-)
-
-test_that("TDM before first dose is considered a true initial value", {
-  expect_equal(sum(is.na(s$y)), 0)
-  expect_equal(s$y[1], 500)
-  expect_equal(round(s$y[3]), 427)
-  expect_equal(round(s$y[13]),1157)
-})
-
-test_that("Variables are set (also in first row) when TDM before first dose", {
-  expect_equal(round(s$CONC[1:5], 1), c(500, 462.2, 427.3, 395.1, 365.3))
-})
diff --git a/tests/testthat/test_time_rounding_bug.R b/tests/testthat/test_time_rounding_bug.R
deleted file mode 100644
index 018af52d..00000000
--- a/tests/testthat/test_time_rounding_bug.R
+++ /dev/null
@@ -1,15 +0,0 @@
-# rounding time should not produce NAs in sim
-## time-rounding bug 20170804
-
-test_that("No NAs related to rounding", {
-  # Uses models defined in setup.R
-  p <- list(CL = 5, V = 50, Q = 10, V2 = 150)
-  r1 <- new_regimen(amt = 100, times = c(0, 24, 36), type = "infusion")
-  dat1 <- sim_ode(
-    ode = mod_2cmt_iv,
-    parameters = p,
-    regimen = r1,
-    t_obs=seq(0, 150, length.out = 100)
-  )
-  expect_equal(sum(is.na(dat1$y)), 0)
-})
diff --git a/tests/testthat/test_timevar_cov.R b/tests/testthat/test_timevar_cov.R
deleted file mode 100644
index b9b9248a..00000000
--- a/tests/testthat/test_timevar_cov.R
+++ /dev/null
@@ -1,111 +0,0 @@
-# Uses model defined in setup.R (conditional, NOT_CRAN only):
-# - mod_2cmt_timevar
-
-skip_on_cran()
-
-par_timevar <- list(CL = 3, V = 50, Q = 2.5, V2 = 70)
-reg_timevar <- new_regimen(
-  amt = 250,
-  n = 60,
-  interval = 6,
-  type = 'infusion',
-  t_inf = 1
-)
-
-test_that("timevarying covariates handled", {
-  # CLi changes by several orders of magnitude after
-  # steady state is achieved, which should produce
-  # a new steady state that is much lower
-  covs <- list(
-    CRCL = new_covariate(
-      value = c(4.662744, 5.798767, 6.195943, 6.2, 600),
-      times = c(0, 18, 29, 209, 210),
-      implementation = "locf"
-    )
-  )
-  t_obs <- seq(0, 360, 0.1)
-  sim1 <- sim_ode(
-    mod_2cmt_timevar,
-    parameters = par_timevar,
-    regimen = reg_timevar,
-    covariates = covs,
-    only_obs = TRUE,
-    t_obs = t_obs,
-    output_include = list(parameters = TRUE, covariates = TRUE)
-  )
-  expect_equal(sim1$CRCL[sim1$t == 209], 6.2)
-  expect_equal(sim1$CRCL[sim1$t == 210], 600)
-  expect_true(sim1$y[sim1$t == 35 * 6] > 10 * sim1$y[sim1$t == 60 * 6])
-})
-
-test_that("timevarying covariates are interpolated and affect PK", {
-  covs_locf <- list(
-    CRCL = new_covariate(
-      times = c(0, 48, 96),
-      value = c(3, 10, 3),
-      implementation = "locf"
-    )
-  )
-  covs_inter <- list(
-    CRCL = new_covariate(
-      times = c(0, 48, 96),
-      value = c(3, 10, 3),
-      implementation = "interpolate"
-    )
-  )
-  t_obs <- seq(0, 120, 2)
-  sim2_inter <- sim_ode(
-    mod_2cmt_timevar,
-    parameters = par_timevar,
-    regimen = reg_timevar,
-    covariates = covs_inter,
-    only_obs = TRUE,
-    t_obs = t_obs,
-    output_include = list(parameters = TRUE, covariates = TRUE, variables = TRUE)
-  )
-  sim2_locf <- sim_ode(
-    mod_2cmt_timevar,
-    parameters = par_timevar,
-    regimen = reg_timevar,
-    covariates = covs_locf,
-    only_obs = TRUE,
-    t_obs = t_obs,
-    output_include = list(parameters = TRUE, covariates = TRUE, variables = TRUE)
-  )
-
-  ## Check covariate changing for inter, but not for locf
-  expect_equal(
-    round(sim2_inter$CRCL, 3)[1:10],
-    c(3, 3.292, 3.583, 3.875, 4.167, 4.458, 4.75, 5.042, 5.333, 5.625)
-  )
-  expect_equal(
-    round(sim2_locf$CRCL, 3)[1:10],
-    rep(3, 10)
-  )
-
-  ## Check interpolated covariates actually affect PK parameters
-  expect_equal(
-    round(sim2_inter$CLi, 3)[1:10],
-    c(6, 6.292, 6.583, 6.875, 7.167, 7.458, 7.75, 8.042, 8.333, 8.625)
-  )
-  expect_equal(
-    round(sim2_locf$CLi, 3)[1:10],
-    rep(6, 10)
-  )
-
-  ## Check interpolated covariates actually affect simulated conc
-  expect_equal(
-    round(sim2_inter$y, 3)[1:10],
-    c(0, 0.32, 0.611, 0.788, 1.205, 1.531, 1.708, 2.098, 2.379, 2.503)
-  )
-  expect_equal(
-    round(sim2_locf$y, 3)[1:10],
-    c(0, 0.321, 0.617, 0.805, 1.239, 1.598, 1.813, 2.251, 2.598, 2.792)
-  )
-
-  ## Visual check:
-  # ggplot(sim2_inter, aes(x = t, y = y)) +
-  #   geom_line() +
-  #   geom_line(data = sim2_locf, colour = "blue")
-
-})

From 76b51487ee06284c074ca1de03d33e897cc51454 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 12:04:41 -0800
Subject: [PATCH 08/10] try fix test

---
 tests/testthat/setup.R    | 7 -------
 tests/testthat/test_sim.R | 8 ++++++--
 2 files changed, 6 insertions(+), 9 deletions(-)

diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
index ed1138f5..1e9878ac 100644
--- a/tests/testthat/setup.R
+++ b/tests/testthat/setup.R
@@ -96,13 +96,6 @@ dose_in_cmt_2 <- new_ode_model(
   cpp_show_code = FALSE
 )
 
-# --- From test_cmt_mapping.R: 1cmt oral with cmt_mapping ---
-pk1cmt_oral_cmt_mapping <- new_ode_model(
-  code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
-  obs = list(cmt = 2, scale = "V"),
-  cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
-)
-
 # --- From test_multi_obs.R: Multi-observation model ---
 vars_multi_obs <- c("CONC", "METAB", "METAB2", "ACT")
 pk_multi_obs <- new_ode_model(
diff --git a/tests/testthat/test_sim.R b/tests/testthat/test_sim.R
index d8984bec..62a29a1e 100644
--- a/tests/testthat/test_sim.R
+++ b/tests/testthat/test_sim.R
@@ -1309,8 +1309,12 @@ describe("Mixture models", {
 })
 
 describe("Compartment mapping", {
-  # Uses model defined in setup.R:
-  # - pk1cmt_oral_cmt_mapping
+
+  pk1cmt_oral_cmt_mapping <- new_ode_model(
+    code = "dAdt[1] = -KA*A[1]; dAdt[2] = KA*A[1] - (CL/V)*A[2];",
+    obs = list(cmt = 2, scale = "V"),
+    cmt_mapping = list(oral = 1, infusion = 2, bolus = 2)
+  )
 
   test_that("Compartment mapping is added to attributes", {
     expect_equal(attr(pk1cmt_oral_cmt_mapping, "cmt_mapping")[["oral"]], 1)

From 8020fe9ad05fa29b2c73e7bba7131ff4bc3e616c Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 16 Jan 2026 12:19:51 -0800
Subject: [PATCH 09/10] add seed

---
 tests/testthat/test_get_var_y.R | 27 +++++++++++++++++----------
 1 file changed, 17 insertions(+), 10 deletions(-)

diff --git a/tests/testthat/test_get_var_y.R b/tests/testthat/test_get_var_y.R
index fb130274..750a5592 100644
--- a/tests/testthat/test_get_var_y.R
+++ b/tests/testthat/test_get_var_y.R
@@ -15,7 +15,7 @@ reg <- new_regimen(
 par <- list(CL = 5, V = 50)
 omega <- c(0.1, 0.0, 0.1)
 t_obs <- c(2, 48)
-res <- sim_ode(
+res <- sim(
   mod_1cmt_iv,
   parameters = par,
   regimen = reg,
@@ -32,7 +32,8 @@ test_that("delta approximation and full simulation match", {
     t_obs = t_obs,
     regimen = reg,
     omega = omega,
-    method = "delta"
+    method = "delta",
+    seed = 12345
   )
   v_sim <- get_var_y(
     model = mod_1cmt_iv,
@@ -42,7 +43,8 @@ test_that("delta approximation and full simulation match", {
     regimen = reg,
     omega = omega,
     method="sim",
-    n_ind = 2500
+    n_ind = 2500,
+    seed = 12345
   )
   expect_true(max(abs(v_delta$regular - v_sim$regular)/v_sim$regular) <0.1)
   expect_true(max(abs(v_delta$log - v_sim$log)/v_sim$log) < 0.15)
@@ -66,7 +68,8 @@ test_that("Confidence interval instead of SD", {
       regimen = reg2,
       omega = omega,
       method = "delta",
-      q = CI_range
+      q = CI_range,
+      seed = 12345
     )
   v2_sim <- get_var_y(
     model = mod_1cmt_iv,
@@ -76,7 +79,8 @@ test_that("Confidence interval instead of SD", {
     omega = omega,
     method ="sim",
     n_ind = 2500,
-    q = CI_range
+    q = CI_range,
+    seed = 12345
   )
   expect_equal(dim(v1_delta$regular), c(2, 2))
   expect_equal(dim(v2_sim$regular), c(2, 2))
@@ -91,7 +95,6 @@ test_that("Confidence interval instead of SD", {
   expect_true(max(abs(v1_delta_v - v2_sim_v )/v2_sim_v) < 0.1)
 })
 
-
 test_that("Two compartment model", {
   skip_on_cran()
   set.seed(80)
@@ -114,7 +117,8 @@ test_that("Two compartment model", {
     parameters = par2,
     t_obs = t_obs,
     regimen = reg,
-    omega = omega2
+    omega = omega2,
+    seed = 12345
   )
   v2 <- get_var_y(
     model = mod_2cmt_iv,
@@ -123,7 +127,8 @@ test_that("Two compartment model", {
     regimen = reg,
     omega = omega2,
     method = "sim",
-    n_ind = 2000
+    n_ind = 2000,
+    seed = 12345
   )
 
   expect_true(all(abs(((v1$regular - v2$regular)/v2$regular)) < 0.05))
@@ -150,7 +155,8 @@ test_that("One compartment with MM kinetics", {
     parameters = par3,
     t_obs = t_obs,
     regimen = reg,
-    omega = omega3
+    omega = omega3,
+    seed = 12345
   )
   v2 <- get_var_y(
     model = mod_1cmt_iv_mm,
@@ -159,7 +165,8 @@ test_that("One compartment with MM kinetics", {
     regimen = reg,
     omega = omega3,
     method = "sim",
-    n_ind = 2000
+    n_ind = 2000,
+    seed = 12345
   )
 
   expect_true(all(abs(v1$regular - v2$regular)/res$y < 0.05))

From 44b839f6f13f738e58dc15e89752f3dc11bc2e5c Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Thu, 22 Jan 2026 14:31:55 -0800
Subject: [PATCH 10/10] add TODO

---
 tests/testthat/test_advan.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/testthat/test_advan.R b/tests/testthat/test_advan.R
index 86435c06..c5bc8b25 100644
--- a/tests/testthat/test_advan.R
+++ b/tests/testthat/test_advan.R
@@ -296,7 +296,7 @@ describe("ADVANs with covariates", {
       times = seq(0, interval * n_days * (24/interval), interval),
       type = "bolus"
     )
-    ## there is slight difference in how bolus doses are handled.
+    ## TODO: 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]