Changed test data generation and single- and multi-individual testing…

… structures for constant model.

tests/testthat/fixtures/constant/make_constant.R

@@ -4,23 +4,12 @@ DE <- function(y, pars){
   return(pars[1])
 }
 
-#Normally distributed error
-add_error_process <- function(y, sigma_e=0.001){
-  return(y + rnorm(length(y), mean=0, sd=sigma_e))
-}
-
 #Function to generate distribution of DE parameters
 DE_par_generation <- function(n_ind, pars=list(mean=0, sd=1)){
   par_sample <- data.frame(beta=exp(rnorm(n_ind, mean = pars[[1]], sd=pars[[2]])))
   return(par_sample)
 }
 
-#Generate initial conditions
-initial_condition_generation <- function(n_ind, pars=list(mean=2, sd=1)){
-  y_0_sample <- data.frame(y_0=exp(rnorm(n_ind, mean = pars[[1]], sd=pars[[2]])))
-  return(y_0_sample)
-}
-
 #Set required values
 model_name <- "constant"
 n_ind <- 3 #Number of individuals for multi-individual data. Single individual takes the first.
@@ -30,14 +19,14 @@ time = seq(from = 0, by = interval, length.out = n_obs_per_ind)
 
 #Produce parameters
 DE_pars <- DE_par_generation(n_ind)
-initial_conditions <- initial_condition_generation(n_ind)
+initial_conditions <- data.frame(y_0=exp(rnorm(n_ind, mean = 2, sd=1)))
 
 #Generate true values
 true_data <- rmot_build_true_test_data(n_ind, n_obs_per_ind, interval,
                                        DE_pars, initial_conditions, DE)
 
 #Generate observations
-y_obs <- add_error_process(true_data$y_true)
+y_obs <- rmot_add_normal_error(true_data$y_true)
 
 #Export datasets
 rmot_export_test_data(n_obs_per_ind,
@@ -48,52 +37,3 @@ rmot_export_test_data(n_obs_per_ind,
                       initial_conditions,
                       true_data,
                       model_name)
-
-#Build model fit single ind
-const_data <- readRDS(test_path("fixtures/constant/constant_data_single_ind.rds"))
-suppressWarnings( #Suppresses stan warnings
-  constant_single_ind_test <- rmot_model("constant_single_ind") |>
-    rmot_assign_data(n_obs = const_data$n_obs, #integer
-                     y_obs = const_data$y_obs,
-                     obs_index = const_data$obs_index, #vector length N_obs
-                     time = const_data$time, #Vector length N_obs
-                     y_0_obs = const_data$y_0_obs #vector length N_ind
-    ) |>
-    rmot_run(chains = 1, iter = 300, verbose = FALSE, show_messages = FALSE, seed=1)
-)
-
-# Extract 100 rows of the samples
-constant_baseline_output_single_ind <- rstan::extract(constant_single_ind_test, permuted = FALSE, inc_warmup = FALSE) |>
-  as.data.frame() |>
-  head(n=100)
-
-# Save the output to compare againest
-saveRDS(constant_baseline_output_single_ind,
-        file = test_path("fixtures/constant/constant_baseline_output_single_ind.rds")
-)
-
-#Build model fit multiple inds
-const_data <- readRDS(test_path("fixtures/constant/constant_data_multi_ind.rds"))
-suppressWarnings( #Suppresses stan warnings
-  constant_multi_ind_test <- rmot_model("constant_multi_ind") |>
-    rmot_assign_data(n_obs = const_data$n_obs, #integer
-                     n_ind = const_data$n_ind, #integer
-                     y_obs = const_data$y_obs, #vector length N_obs
-                     obs_index = const_data$obs_index, #vector length N_obs
-                     time = const_data$time, #Vector length N_obs
-                     ind_id = const_data$ind_id, #Vector length N_obs
-                     y_0_obs = const_data$y_0_obs #vector length N_ind
-    ) |>
-    rmot_run(chains = 1, iter = 300, verbose = FALSE, show_messages = FALSE, seed=1)
-)
-
-# Extract 100 rows of the samples
-constant_baseline_output_multi_ind <- rstan::extract(constant_multi_ind_test, permuted = FALSE, inc_warmup = FALSE) |>
-  as.data.frame() |>
-  head(n=100)
-
-# Save the output to compare againest
-saveRDS(constant_baseline_output_multi_ind,
-        file = test_path("fixtures/constant/constant_baseline_output_multi_ind.rds")
-        )
-

tests/testthat/helper.R → tests/testthat/helper-data_generation.R

@@ -1,6 +1,4 @@
-
 #General functions for data generation
-
 #Runge-Kutta 4th order
 rmot_rk4_est <- function(y_0, DE, pars, step_size, n_step){
   runge_kutta_int <- c(y_0)
@@ -46,6 +44,10 @@ rmot_build_true_test_data <- function(n_ind, n_obs, interval,
   return(true_data)
 }
 
+rmot_add_normal_error <- function(y, sigma_e=0.001){
+  return(y + rnorm(length(y), mean=0, sd=sigma_e))
+}
+
 #Save data to files
 rmot_export_test_data <- function(n_obs_per_ind,
                                   n_ind,
@@ -63,6 +65,7 @@ rmot_export_test_data <- function(n_obs_per_ind,
     obs_index = 1:n_obs_per_ind, #Vector indexed by n_obs
     time = time, #Vector indexed by n_obs
     y_0_obs = y_obs[1], #Number
+    n_pars = ncol(DE_pars), #Number
     single_true_data = list(
       DE_pars = DE_pars[1,],
       initial_conditions = initial_conditions[1,],
@@ -80,6 +83,7 @@ rmot_export_test_data <- function(n_obs_per_ind,
     time = rep(time, times=n_ind), #Vector indexed by n_obs
     ind_id = sort(rep(1:n_ind, times = n_obs_per_ind)), #Vector indexed by n_obs
     y_0_obs = y_obs[seq(from = 1, to=n_ind*n_obs_per_ind, by=n_obs_per_ind)], #Vector indexed by n_ind
+    n_pars = ncol(DE_pars),
     multi_true_data = list(
       DE_pars = DE_pars,
       initial_conditions = initial_conditions,
@@ -88,7 +92,7 @@ rmot_export_test_data <- function(n_obs_per_ind,
     )
   )
 
-  if(! exists(test_path("fixtures", model_name))) dir.create(test_path("fixtures", model_name))
+  if(! dir.exists(test_path("fixtures", model_name))) dir.create(test_path("fixtures", model_name))
   filename <- paste("tests/testthat/fixtures", "/", model_name, "/", model_name, "_data", sep="")
   saveRDS(single_ind_data, file=paste(filename, "single_ind.rds", sep="_"))
   saveRDS(multi_ind_data, file=paste(filename, "multi_ind.rds", sep="_"))

tests/testthat/helper-testing.R

@@ -0,0 +1,72 @@
+rmot_test_model_functions <- function(model_name){
+  single_ind_name <- paste0(model_name, "_single_ind")
+  multi_ind_name <- paste0(model_name, "_multi_ind")
+
+  # Single individual
+  expect_named(rmot_model(single_ind_name))
+  expect_type(rmot_model(single_ind_name), "list")
+  expect_visible(rmot_model(single_ind_name))
+  #Multiple individuals
+  expect_named(rmot_model(multi_ind_name))
+  expect_type(rmot_model(multi_ind_name), "list")
+  expect_visible(rmot_model(multi_ind_name))
+}
+
+rmot_test_single_individual <- function(model_name,
+                                        par_names){
+  data <- readRDS(test_path("fixtures", model_name,
+                            paste0(model_name, "_data_single_ind.rds")))
+
+  # Test constant single individual
+  suppressWarnings( #Suppresses stan warnings
+    single_ind_test <- rmot_model(paste0(model_name, "_single_ind")) |>
+      rmot_assign_data(n_obs = data$n_obs, #integer
+                       y_obs = data$y_obs,
+                       obs_index = data$obs_index, #vector length N_obs
+                       time = data$time, #Vector length N_obs
+                       y_0_obs = data$y_0_obs #vector length N_ind
+      ) |>
+      rmot_run(chains = 1, iter = 1000, verbose = FALSE, show_messages = FALSE)
+  )
+
+  # Extract samples and check if parameter estimates are reasonable.
+  ind_samples <- rstan::extract(single_ind_test, permuted = TRUE,
+                                         inc_warmup = FALSE)
+  par_ests <- c()
+  for(i in length(par_names)){
+    par_ests[i] <- mean(ind_samples[[par_names[i]]])
+  }
+
+  initial_condition <- mean(ind_samples$ind_y_0)
+  expect_equal(par_ests, data$single_true_data$DE_pars,
+               tolerance = 1e-1)
+  expect_equal(initial_condition, data$single_true_data$initial_conditions,
+               tolerance = 1e-1)
+
+  # hecks for output existence and type
+  expect_visible(single_ind_test)
+  expect_s4_class(single_ind_test, "stanfit")
+}
+
+rmot_test_multi_individual <- function(model_name, data, est_dim){
+  # Test constant multi-individual
+  suppressWarnings( #Suppresses stan warnings
+    multi_ind_test <- rmot_model(paste0(model_name, "_multi_ind")) |>
+      rmot_assign_data(n_obs = data$n_obs, #integer
+                       n_ind = data$n_ind, #integer
+                       y_obs = data$y_obs, #vector length N_obs
+                       obs_index = data$obs_index, #vector length N_obs
+                       time = data$time, #Vector length N_obs
+                       ind_id = data$ind_id, #Vector length N_obs
+                       y_0_obs = data$y_0_obs #vector length N_ind
+      ) |>
+      rmot_run(chains = 2, iter = 100, verbose = FALSE, show_messages = FALSE)
+  )
+
+  # Extract samples
+  multi_samples <- rstan::extract(multi_ind_test, permuted = FALSE, inc_warmup = TRUE)
+  expect_equal(dim(multi_samples), c(100, 2, est_dim))
+
+  expect_visible(multi_ind_test)
+  expect_s4_class(multi_ind_test, "stanfit")
+}

tests/testthat/test-rmot_models.R

@@ -1,17 +1,7 @@
-test_that("Model structures", {
+test_that("Model structures: Linear", {
   expect_named(rmot_model("linear"))
   expect_type(rmot_model("linear"), "list")
   expect_visible(rmot_model("linear"))
-
-  # Constant models
-  # Single individual
-  expect_named(rmot_model("constant_single_ind"))
-  expect_type(rmot_model("constant_single_ind"), "list")
-  expect_visible(rmot_model("constant_single_ind"))
-  #Multiple individuals
-  expect_named(rmot_model("constant_multi_ind"))
-  expect_type(rmot_model("constant_multi_ind"), "list")
-  expect_visible(rmot_model("constant_multi_ind"))
 })
 
 test_that("Execution and output: Linear", {
@@ -37,66 +27,3 @@ test_that("Execution and output: Linear", {
   expect_visible(lm_test)
   expect_s4_class(lm_test, "stanfit")
 })
-
-test_that("Execution: Constant single individual", {
-  const_data <- readRDS(test_path("fixtures", "constant",
-                                   "constant_data_single_ind.rds"))
-  const_single_ind_baseline_output <- readRDS(test_path("fixtures", "constant",
-                                                        "constant_baseline_output_single_ind.rds"))
-
-   # Test constant single individual
-  suppressWarnings( #Suppresses stan warnings
-    constant_single_ind_test <- rmot_model("constant_single_ind") |>
-      rmot_assign_data(n_obs = const_data$n_obs, #integer
-                       y_obs = const_data$y_obs,
-                       obs_index = const_data$obs_index, #vector length N_obs
-                       time = const_data$time, #Vector length N_obs
-                       y_0_obs = const_data$y_0_obs #vector length N_ind
-      ) |>
-      rmot_run(chains = 1, iter = 300, verbose = FALSE, show_messages = FALSE,
-      seed = 1)
-  )
-
-  # Extract samples
-  constant_ind_samples <- rstan::extract(constant_single_ind_test, permuted = FALSE, inc_warmup = FALSE) |>
-    as.data.frame() |>
-    head(n=100)
-
-  expect_equal(constant_ind_samples,
-               const_single_ind_baseline_output, tolerance = 1e-5)
-  expect_visible(constant_single_ind_test)
-  expect_s4_class(constant_single_ind_test, "stanfit")
-
-})
-
-test_that("Execution: Constant multiple individuals", {
-  const_data <- readRDS(test_path("fixtures", "constant",
-                                  "constant_data_multi_ind.rds"))
-  const_multi_ind_baseline_output <- readRDS(test_path("fixtures", "constant",
-                                                       "constant_baseline_output_multi_ind.rds"))
-
-  # Test constant multi-individual
-  suppressWarnings( #Suppresses stan warnings
-    constant_multi_ind_test <- rmot_model("constant_multi_ind") |>
-      rmot_assign_data(n_obs = const_data$n_obs, #integer
-                       n_ind = const_data$n_ind, #integer
-                       y_obs = const_data$y_obs, #vector length N_obs
-                       obs_index = const_data$obs_index, #vector length N_obs
-                       time = const_data$time, #Vector length N_obs
-                       ind_id = const_data$ind_id, #Vector length N_obs
-                       y_0_obs = const_data$y_0_obs #vector length N_ind
-      ) |>
-      rmot_run(chains = 1, iter = 300, verbose = FALSE, show_messages = FALSE,
-      seed = 1)
-  )
-
-  # Extract samples
-  constant_multi_samples <- rstan::extract(constant_multi_ind_test, permuted = FALSE, inc_warmup = FALSE) |>
-    as.data.frame() |>
-    head(n=100)
-
-  expect_equal(constant_multi_samples,
-               const_multi_ind_baseline_output, tolerance = 1e-5)
-  expect_visible(constant_multi_ind_test)
-  expect_s4_class(constant_multi_ind_test, "stanfit")
-})

tests/testthat/test-rmot_models_constant.R

@@ -0,0 +1,36 @@
+#Testing for constant model
+test_that("Model structures: Constant", {
+  # Single individual
+  expect_named(rmot_model("constant_single_ind"))
+  expect_type(rmot_model("constant_single_ind"), "list")
+  expect_visible(rmot_model("constant_single_ind"))
+  #Multiple individuals
+  expect_named(rmot_model("constant_multi_ind"))
+  expect_type(rmot_model("constant_multi_ind"), "list")
+  expect_visible(rmot_model("constant_multi_ind"))
+})
+
+test_that("Execution: Constant single individual", {
+  model_name <- "constant"
+  par_names <- "ind_beta"
+
+  rmot_test_single_individual(model_name, par_names)
+})
+
+test_that("Execution: Constant multiple individuals", {
+  model_name <- "constant"
+
+  data <- readRDS(test_path("fixtures", "constant",
+                            "constant_data_multi_ind.rds"))
+
+  #Dimension is:
+  est_dim <- data$n_ind +     #Initial condition
+    data$n_pars * data$n_ind + #Individual parameters
+    data$n_pars * 2  +         #Population parameters
+    1 +                        #Global error
+    data$n_obs +               #y_ij
+    data$n_obs +               #Delta y_ij
+    1                          #lp__
+
+  rmot_test_multi_individual(model_name, data, est_dim)
+})