Dragons rework

DESCRIPTION

@@ -43,7 +43,9 @@ Suggests:
     here,
     patchwork,
     deSolve,
-    cowplot
+    cowplot,
+    mixtools,
+    MASS
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 LazyData: true

NAMESPACE

@@ -1,10 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
+export(hmde_affine_de)
 export(hmde_assign_data)
 export(hmde_canham_de)
 export(hmde_const_de)
 export(hmde_extract_estimates)
-export(hmde_linear_de)
 export(hmde_model)
 export(hmde_model_des)
 export(hmde_model_names)

R/hmde_assign_data.R

@@ -46,7 +46,7 @@ hmde_assign_data <- function(model_template, data = NULL,...){
   for(i in model_fields){ # Iterate through required fields and fill them
     if(i %in% user_fields){
       model_template <- purrr::list_modify(model_template, !!!data[i])
-    } else {
+    } else if(is.null(model_template[[i]])){ # allows for default values
       model_template[[i]] <- switch(
         i,
         n_obs = length(data$y_obs),

R/hmde_model_des.R

@@ -17,7 +17,7 @@ hmde_model_des <- function(model = NULL){
     canham_multi_ind = hmde_canham_de,
     vb_single_ind = hmde_vb_de,
     vb_multi_ind = hmde_vb_de,
-    linear_single_ind = hmde_linear_de
+    affine_single_ind = hmde_affine_de
   )
 
   return(output)
@@ -65,15 +65,15 @@ hmde_vb_de <- function(y = NULL, pars = NULL){
     )
 }
 
-#' Differential equation for linear growth single individual model
+#' Differential equation for affine growth single individual model
 #' @param y input real
 #' @param pars list of parameters beta_0, beta_1
 #'
 #' @return value of differential equation at y
 #' @export
 
-hmde_linear_de <- function(y = NULL, pars = NULL){
+hmde_affine_de <- function(y = NULL, pars = NULL){
   return(
-    pars[[1]] + pars[[2]] * y
+    pars[[1]] - pars[[2]] * y
     )
 }

R/hmde_model_names.R

@@ -10,7 +10,7 @@ hmde_model_names <- function(){
               "canham_multi_ind",
               "vb_single_ind",
               "vb_multi_ind",
-              "linear_single_ind")
+              "affine_single_ind")
 
   return(output)
 }

R/hmde_model_pars.R

@@ -18,7 +18,7 @@ hmde_model_pars <- function(model=NULL){
                    canham_multi_ind = hmde_canham_multi_ind_pars(),
                    vb_single_ind = hmde_vb_single_ind_pars(),
                    vb_multi_ind = hmde_vb_multi_ind_pars(),
-                   linear_single_ind = hmde_linear_single_ind_pars())
+                   affine_single_ind = hmde_affine_single_ind_pars())
 
   return(output)
 }
@@ -95,13 +95,13 @@ hmde_vb_multi_ind_pars <- function(){
        model = "vb_multi_ind")
 }
 
-#' Parameter names for linear growth single individual model
+#' Parameter names for affine growth single individual model
 #' @keywords internal
 #' @noRd
 #'
-hmde_linear_single_ind_pars <- function(){
+hmde_affine_single_ind_pars <- function(){
   list(measurement_pars_names = c("y_hat"),
        individual_pars_names = c("ind_beta_0", "ind_beta_1"),
        error_pars_names = c(NULL),
-       model = "linear_single_ind")
+       model = "affine_single_ind")
 }

R/hmde_models.R

@@ -18,7 +18,7 @@ hmde_model <- function(model=NULL){
                    canham_multi_ind = hmde_canham_multi_ind(),
                    vb_single_ind = hmde_vb_single_ind(),
                    vb_multi_ind = hmde_vb_multi_ind(),
-                   linear_single_ind = hmde_linear_single_ind())
+                   affine_single_ind = hmde_affine_single_ind())
 
   class(output) <- "hmde_object"
 
@@ -106,16 +106,19 @@ hmde_vb_multi_ind <- function(){
        model = "vb_multi_ind")
 }
 
-#' Data configuration template for linear growth single individual model
+#' Data configuration template for affine growth single individual model
 #' @keywords internal
 #' @noRd
 #'
-hmde_linear_single_ind <- function(){
+hmde_affine_single_ind <- function(){
   list(step_size = NULL,
        n_obs = NULL,
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
        y_bar = NULL,
-       model = "linear_single_ind")
+       int_method = NULL,
+       prior_means = c(1,1),
+       prior_sds = c(2,2),
+       model = "affine_single_ind")
 }

R/hmde_run.R

@@ -20,7 +20,7 @@ hmde_run <- function(model_template, ...) {
          canham_multi_ind = rstan::sampling(stanmodels$canham_multi_ind, data = model_template, ...),
          vb_single_ind = rstan::sampling(stanmodels$vb_single_ind, data = model_template, ...),
          vb_multi_ind = rstan::sampling(stanmodels$vb_multi_ind, data = model_template, ...),
-         linear_single_ind = rstan::sampling(stanmodels$linear_single_ind, data = model_template, ...))
+         affine_single_ind = rstan::sampling(stanmodels$affine_single_ind, data = model_template, ...))
 
   return(out)
 }

R/stanmodels.R

@@ -1,14 +1,14 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "linear_single_ind", "vb_multi_ind", "vb_single_ind")
+stanmodels <- c("affine_single_ind", "canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "vb_multi_ind", "vb_single_ind")
 
 # load each stan module
+Rcpp::loadModule("stan_fit4affine_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4canham_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4canham_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4constant_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4constant_single_ind_mod", what = TRUE)
-Rcpp::loadModule("stan_fit4linear_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4vb_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4vb_single_ind_mod", what = TRUE)
 

hmde.Rproj

@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: 276917f4-3da1-492e-bb45-7ffd205ea048
 
 RestoreWorkspace: No
 SaveWorkspace: No

inst/stan/linear_single_ind.stan → inst/stan/affine_single_ind.stan

@@ -45,6 +45,16 @@ functions{
 
     return y_hat;
   }
+
+  vector DE_RK45(real t, vector y, real ind_const, real beta_1, real y_bar){
+    vector[size(y)] dydt = ind_const - (beta_1 * (y-y_bar));
+    return dydt;
+  }
+
+  real analytic_solution(real t, real y_0, real ind_const, real beta_1, real y_bar){
+    real y_0_translate = y_0 - y_bar;
+    return ind_const/beta_1 + (y_0_translate - (ind_const/beta_1)) * exp(-beta_1 * t) + y_bar;
+  }
 }
 
 // Data structure
@@ -55,6 +65,9 @@ data {
   int obs_index[n_obs];
   real time[n_obs];
   real y_bar;
+  int int_method;
+  real prior_means[2]; #vector of means for beta parameter priors
+  real prior_sds[2]; #Vector of SDs for beta parameter priors
 }
 
 // The parameters accepted by the model.
@@ -72,6 +85,7 @@ parameters {
 model {
   real y_hat[n_obs];
   array[3] real pars;
+  vector[1] y_temp;
 
   pars[1] = ind_const;
   pars[2] = ind_beta_1;
@@ -85,7 +99,24 @@ model {
 
     if(i < n_obs){
       //Estimate next size
-      y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      if(int_method == 1){ //RK4
+        y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      }
+
+      if(int_method == 2){ //RK45
+        y_temp[1] = y_hat[i];
+        y_hat[i+1] = ode_rk45(DE_RK45, y_temp,
+                      time[i], {time[i+1]},
+                      ind_const, ind_beta_1, y_bar)[1][1];
+      }
+
+      if(int_method == 3){ //Analytic solution
+        y_hat[i+1] = analytic_solution((time[i+1] - time[1]),
+                          ind_y_0,
+                          ind_const,
+                          ind_beta_1,
+                          y_bar);
+      }
     }
   }
 
@@ -94,14 +125,16 @@ model {
 
   //Priors
   //Individual level
-  ind_const ~lognormal(0, 2);
-  ind_beta_1 ~lognormal(0, 2);
+  ind_const ~lognormal(log(prior_means[1]), prior_sds[1]);
+  ind_beta_1 ~lognormal(log(prior_means[2]), prior_sds[2]);
 }
 
 generated quantities{
   real y_hat[n_obs];
   array[3] real pars;
   real ind_beta_0;
+  vector[1] y_temp;
+  int vers = 1;
 
   ind_beta_0 = ind_const + ind_beta_1*y_bar;
 
@@ -117,7 +150,24 @@ generated quantities{
 
     if(i < n_obs){
       //Estimate next size
-      y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      if(int_method == 1){ //RK4
+        y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      }
+
+      if(int_method == 2){ //RK45
+        y_temp[1] = y_hat[i];
+        y_hat[i+1] = ode_rk45(DE_RK45, y_temp,
+                      time[i], {time[i+1]},
+                      ind_const, ind_beta_1, y_bar)[1][1];
+      }
+
+      if(int_method == 3){ //Analytic solution
+        y_hat[i+1] = analytic_solution((time[i+1] - time[1]),
+                          ind_y_0,
+                          ind_const,
+                          ind_beta_1,
+                          y_bar);
+      }
     }
   }
 }

tests/testthat/helper-data_generation.R

@@ -60,7 +60,7 @@ hmde_export_test_data <- function(n_obs_per_ind,
                                   true_data,
                                   model_name,
                                   step_size = 1){
-  if(grepl("linear", model_name)){ #Step size for linear model
+  if(grepl("affine", model_name)){ #Step size for affine model
     single_ind_data <- list(
       step_size = step_size, #Number for model RK4 alg
       n_obs = n_obs_per_ind, #Number

tests/testthat/test-hmde_assign_data.R

@@ -119,3 +119,33 @@ test_that("Execution and output: bad input", {
       hmde_assign_data(data = mtcars)
   )
 })
+
+
+test_that("Execution and output: default values", {
+  time <- 1:5
+  y_obs <- 1:5
+  obs_index <- 1:5
+
+  default_used <- hmde_model("affine_single_ind") |>
+    hmde_assign_data(n_obs = length(y_obs),
+                     y_obs= y_obs,
+                     obs_index = obs_index,
+                     time = time,
+                     y_bar = mean(y_obs),
+                     step_size = 1,
+                     int_method = 1)
+
+  expect_equal(default_used$prior_means, c(1,1))
+
+  value_supplied <- hmde_model("affine_single_ind") |>
+    hmde_assign_data(n_obs = length(y_obs),
+                     y_obs= y_obs,
+                     obs_index = obs_index,
+                     time = time,
+                     y_bar = mean(y_obs),
+                     step_size = 1,
+                     int_method = 1,
+                     prior_means = c(5,5))
+
+  expect_equal(value_supplied$prior_means, c(5,5))
+})

tests/testthat/test-hmde_models_affine.R

@@ -0,0 +1,10 @@
+#Testing for affine model
+test_that("Model structures: affine", {
+  # Single individual
+  single_model <- hmde_model("affine_single_ind")
+  expect_named(single_model, c("step_size", "n_obs", "y_obs", "obs_index",
+                               "time", "y_bar", "int_method", "prior_means",
+                               "prior_sds", "model"))
+  expect_type(single_model, "list")
+  expect_visible(single_model)
+})