RXR-2394: support estimation of lagtime

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: PKPDmap
 Type: Package
 Title: MAP Bayesian estimates
-Version: 1.0.3
-Date: 2025-06-18
+Version: 1.1.4
+Date: 2025-07-12
 Author: Ron Keizer, Jasmine Hughes, Kara Woo
 Maintainer: Ron Keizer <ron@insight-rx.com>
 Description: Obtain MAP Bayesian estimates based on individual data and

R/calc_residuals.R

@@ -18,6 +18,7 @@ calc_residuals <- function(
   parameters_population,
   covariates,
   regimen,
+  lagtime,
   omega_full,
   error,
   weights,
@@ -60,6 +61,7 @@ calc_residuals <- function(
       iov_bins = iov_bins,
       output_include = output_include,
       t_init = t_init,
+      lagtime = lagtime,
       ...
     )
   })
@@ -79,6 +81,7 @@ calc_residuals <- function(
       iov_bins = iov_bins,
       A_init = A_init_population,
       t_init = t_init,
+      lagtime = lagtime,
       ...
     )
   })

R/get_individual_parameters_from_fit.R

@@ -1,7 +1,8 @@
 #' Get the individual parameters from the fitted coefficients (etas)
 #' 
+#' @inheritParams get_map_estimates
 #' @param fit fit object
-#' @param parameters list of population parameters
+#' @param nonfixed vector of parameters that are not fixed (i.e. estimated)
 #'
 get_individual_parameters_from_fit <- function(
   fit,

R/get_map_estimates.R

@@ -18,6 +18,8 @@
 #' PKPDsim.
 #' @param error residual error, specified as list with arguments `add` and/or 
 #' `prop` specifying the additive and proportional parts
+#' @param lagtime vector of lagtimes for each compartment in the model, either
+#' numeric vector, or vector of parameters.
 #' @param ltbs log-transform both sides? (`NULL` by default, meaning that it 
 #' will be picked up from the PKPDsim model. Can be overridden with `TRUE`). 
 #' Note: `error` should commonly only have additive part.
@@ -39,7 +41,12 @@
 #' @param ll_func likelihood function, default is `ll_func_PKPDsim` as included 
 #' in this package.
 #' @param optimizer optimization library to use, default is `optim`
-#' @param method optimization method, default `BFGS`
+#' @param method optimization method, default `BFGS`. This method is a
+#' gradient-based method, in which the gradients are computed using finite-
+#' difference method. When the model contains a step-function where a
+#' estimated parameter is involved in the step-function (such as in the case of
+#' lagtime), it is advised to use a non-gradient-based estimation method
+#' such as `Nelder-Mead` to avoid estimation failures.
 #' @param control list of options passed to `optim()` function
 #' @param allow_obs_before_dose allow observation before first dose?
 #' @param type estimation type, options are `map`, `ls`, and `np_hybrid`
@@ -112,6 +119,7 @@ get_map_estimates <- function(
                       include_omega = TRUE,
                       include_error = TRUE,
                       regimen = NULL,
+                      lagtime = NULL,
                       t_init = 0,
                       int_step_size = 0.1,
                       ll_func = ll_func_PKPDsim,
@@ -137,6 +145,9 @@ get_map_estimates <- function(
   if(weight_prior_var == 0) {
     calc_ofv <- calc_ofv_ls
   }
+
+  ## Make sure lagtime is properly formatted (potentially pulled from model)
+  lagtime <- PKPDsim:::parse_lagtime(lagtime, model, parameters)
 
   ## Misc checks:
   check_inputs(model, data, parameters, omega, regimen, censoring, type)
@@ -211,6 +222,7 @@ get_map_estimates <- function(
       n_ind = 1,
       int_step_size = int_step_size,
       regimen = regimen,
+      lagtime = lagtime,
       t_obs = t_obs,
       obs_type = data$obs_type,
       checks = FALSE,
@@ -250,6 +262,7 @@ get_map_estimates <- function(
           t_obs = t_obs,
           model = model,
           regimen = regimen,
+          lagtime = lagtime,
           error = error,
           omega_full = omega$est / weight_prior_var,
           omega_inv = solve(omega$est / weight_prior_var),
@@ -305,6 +318,7 @@ get_map_estimates <- function(
           t_obs = t_obs,
           model = model,
           regimen = regimen,
+          lagtime = lagtime,
           error = error,
           nonfixed = omega$nonfixed,
           transf = transf,
@@ -353,6 +367,7 @@ get_map_estimates <- function(
       obj = obj,
       model = model,
       regimen = regimen,
+      lagtime = lagtime,
       data = data,
       covariates = covariates,
       weights = weights,
@@ -399,6 +414,7 @@ get_map_estimates <- function(
       covariates = covariates,
       int_step_size = int_step_size,
       regimen = regimen,
+      lagtime = lagtime,
       omega_full = omega$full,
       error = error,
       weights = weights,

R/ll_func_PKPDsim.R

@@ -1,29 +1,19 @@
 #' Likelihood function for MAP optimization using PKPDsim
 #' 
-#' @param data vector of
+#' @inheritParams get_map_estimates
 #' @param sim_object design (event-table) obtained from PKPDsim to be used in simulations/optimizations
-#' @param parameters parameter list
-#' @param covariates covariates list
 #' @param nonfixed non-fixed (i.e. estimated) parameters
-#' @param error error model to use, e.g. `list(add = .5, prop = .15)`
-#' @param model PKPDsim model
 #' @param omega_full full omega matrix
 #' @param omega_inv inverse of omega matrix. Passed to this function to avoid doing computations in each iteration of the search.
 #' @param omega_eigen eigenvalue decomposation (logged) of omega matrix. Passed to this function to avoid doing computations in each iteration of the search.
 #' @param sig signficance, used in optimization
-#' @param weights weights for data, generally a vector of weights between 0 and 1.
 #' @param transf transformation function for data, if needed. E.g. `log(x)`. Default is `function(x) = x`.
-#' @param as_eta implement as regular eta instead of exponential eta, can be vector.
 #' @param censoring_idx censoring indices, used for <LOQ data.
 #' @param censoring_label censoring label, used for <LOQ data
-#' @param t_init init time for simulations, default 0.
-#' @param iov_bins IOV bins object
 #' @param calc_ofv function to calculate OFV based on simulated data and set of parameters and omega matrix
-#' @param regimen PKPDsim regimen
 #' @param steady_state_analytic list object with settings for steady state MAP estimation.
 #' @param include_omega include omega in calculation of OFV?
 #' @param include_error include residual error in calculation of OFV?
-#' @param verbose verbose output?
 #' @param eta01 eta1 
 #' @param eta02 eta2 
 #' @param eta03 eta3 
@@ -75,6 +65,7 @@ ll_func_PKPDsim <- function(
   t_init = 0,
   calc_ofv,
   regimen,
+  lagtime = c(0),
   steady_state_analytic,
   include_omega,
   include_error,
@@ -109,7 +100,9 @@ ll_func_PKPDsim <- function(
     sim_object,
     ode = model,
     duplicate_t_obs = TRUE,
-    t_init = t_init)$y)
+    t_init = t_init,
+    lagtime = lagtime
+  )$y)
   dv <- transf(data$y)
   obs_type <- data$obs_type
   ofv_cens <- NULL

R/ll_func_generic.R

@@ -2,6 +2,34 @@
 #' 
 #' @inheritParams ll_func_PKPDsim
 #'
+#' @param omega_full full omega matrix
+#' @param sig signficance, used in optimization
+#' @param eta01 eta1 
+#' @param eta02 eta2 
+#' @param eta03 eta3 
+#' @param eta04 eta4 
+#' @param eta05 eta5 
+#' @param eta06 eta6 
+#' @param eta07 eta7 
+#' @param eta08 eta8 
+#' @param eta09 eta9 
+#' @param eta10 eta10 
+#' @param eta11 eta11 
+#' @param eta12 eta12 
+#' @param eta13 eta13 
+#' @param eta14 eta14 
+#' @param eta15 eta15 
+#' @param eta16 eta16 
+#' @param eta17 eta17 
+#' @param eta18 eta19
+#' @param eta19 eta19
+#' @param eta20 eta20
+#' @param eta21 eta21
+#' @param eta22 eta22
+#' @param eta23 eta23
+#' @param eta24 eta24 
+#' @param ... passed on to PKPDsim
+#' 
 ll_func_generic <- function(
   data,
   # unfortunately seems no other way to do this...
@@ -11,6 +39,7 @@ ll_func_generic <- function(
   parameters,
   covariates = NULL,
   regimen = regimen,
+  lagtime = NULL,
   omega_full = omega_full,
   error = error,
   model,

R/parse_input_data.R

@@ -7,6 +7,12 @@ parse_input_data <- function(
   cols = list(x = "t", y = "y"),
   obs_type_label = NULL
 ) {
+  ## Parse data to include label for obs_type
+  if(is.null(obs_type_label)) {
+    data$obs_type <- 1
+  } else {
+    data$obs_type <- data[[obs_type_label]]
+  }
   colnames(data) <- tolower(colnames(data))
   if(!all(tolower(unlist(cols)) %in% names(data))) {
     stop("Expected column names were not found in data. Please use 'cols' argument to specify column names for independent and dependent variable.")
@@ -18,12 +24,6 @@ parse_input_data <- function(
       data$evid <- 0
     }
   }
-  ## Parse data to include label for obs_type
-  if(is.null(obs_type_label)) {
-    data$obs_type <- 1
-  } else {
-    data$obs_type <- data[[obs_type_label]]
-  }
   if("evid" %in% colnames(data)) {
     data <- data[data$evid == 0,]
   }

R/parse_weights.R

@@ -1,6 +1,7 @@
 #' Parse weights argument
 #' 
 #' @inheritParams get_map_estimates
+#' @param t_obs vector of observation times
 #' 
 parse_weights <- function(weights, t_obs) {
   if(!is.null(weights)) {