From f3018dca1281021d60eb16532ebb0d6eabd33e4b Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Mon, 7 Jul 2025 07:23:27 -0700 Subject: [PATCH 1/6] initial commit --- R/calc_residuals.R | 3 ++ R/get_map_estimates.R | 9 ++++ R/ll_func_PKPDsim.R | 5 ++- R/ll_func_generic.R | 1 + tests/testthat/setup.R | 10 +++++ .../testthat/test-get_map_estimates.lagtime.R | 44 +++++++++++++++++++ 6 files changed, 71 insertions(+), 1 deletion(-) create mode 100644 tests/testthat/setup.R create mode 100644 tests/testthat/test-get_map_estimates.lagtime.R diff --git a/R/calc_residuals.R b/R/calc_residuals.R index 46e3b19..b3c4173 100644 --- a/R/calc_residuals.R +++ b/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, ... ) }) diff --git a/R/get_map_estimates.R b/R/get_map_estimates.R index 046e833..3c44605 100755 --- a/R/get_map_estimates.R +++ b/R/get_map_estimates.R @@ -112,6 +112,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 +138,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) ## Misc checks: check_inputs(model, data, parameters, omega, regimen, censoring, type) @@ -211,6 +215,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 +255,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 +311,7 @@ get_map_estimates <- function( t_obs = t_obs, model = model, regimen = regimen, + lagtime = lagtime, error = error, nonfixed = omega$nonfixed, transf = transf, @@ -353,6 +360,7 @@ get_map_estimates <- function( obj = obj, model = model, regimen = regimen, + lagtime = lagtime, data = data, covariates = covariates, weights = weights, @@ -399,6 +407,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, diff --git a/R/ll_func_PKPDsim.R b/R/ll_func_PKPDsim.R index 8028aa1..e595bf0 100644 --- a/R/ll_func_PKPDsim.R +++ b/R/ll_func_PKPDsim.R @@ -75,6 +75,7 @@ ll_func_PKPDsim <- function( t_init = 0, calc_ofv, regimen, + lagtime = c(0), steady_state_analytic, include_omega, include_error, @@ -109,7 +110,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 diff --git a/R/ll_func_generic.R b/R/ll_func_generic.R index e739a82..e8c42b8 100755 --- a/R/ll_func_generic.R +++ b/R/ll_func_generic.R @@ -11,6 +11,7 @@ ll_func_generic <- function( parameters, covariates = NULL, regimen = regimen, + lagtime = NULL, omega_full = omega_full, error = error, model, diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R new file mode 100644 index 0000000..3fbef06 --- /dev/null +++ b/tests/testthat/setup.R @@ -0,0 +1,10 @@ +mod_1cmt_oral_lagtime <- PKPDsim::new_ode_model( + code = " + dAdt[0] = -KA * A[0] + dAdt[1] = +KA * A[0] -(CL/V) * A[1] + ", + lagtime = c("TLAG", 0), + obs = list(cmt = 2, scale = "V"), + dose = list(cmt = 1, bioav = 1), + parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83) +) \ No newline at end of file diff --git a/tests/testthat/test-get_map_estimates.lagtime.R b/tests/testthat/test-get_map_estimates.lagtime.R new file mode 100644 index 0000000..3bf5c34 --- /dev/null +++ b/tests/testthat/test-get_map_estimates.lagtime.R @@ -0,0 +1,44 @@ +test_that("MAP fit works with oral model with lagtime", { + ## Simulate some data + reg <- PKPDsim::new_regimen( + amt = 100, + times = c(0, 24), + type = "bolus" + ) + data <- PKPDsim::sim( + mod_1cmt_oral_lagtime, + regimen = reg, + parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83), + t_obs = c(0.25, 0.5, 0.75, 1.0, 2, 5, 8, 16), + only_obs = TRUE + ) + ## check that fit works for vector of lagtimes + lagtimes <- c(0.25, 0.5, 0.75, 1.0, 2.0) + est <- c() + for(i in seq(lagtimes)) { + fit <- get_map_estimates( + model = mod_1cmt_oral_lagtime, + parameters = list(CL = 6, V = 40, KA = 0.6, TLAG = lagtimes[i]), # slightly tweaked + regimen = reg, + omega = c( + 0.1, + 0, 0.1, + 0, 0, 0.1 + ), + fixed = c("KA"), + error = list(prop = 0.1, add = 0.15), + data = data, + ## !important! due to lagtime discontinuity, cannot use gradient-based optimizers such as + ## BGFS or L-BFGS-B, eventhough gradients are based on difference method. Optimization + ## with these methods sometimes still work if lagtime initial estimate is higher than + ## true value, but better to use a gradient-free method such as Nelder-Mead. + method = "Nelder-Mead", + verbose = F + ) + est <- c(est, fit$parameters$TLAG) + } + expect_equal( + round(est, 3), + c(0.594, 0.757, 0.856, 1.043, 1.307) + ) +}) From 4af5715b10c37487997c1ba19c8e2d3268951b7a Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Tue, 8 Jul 2025 01:01:18 -0700 Subject: [PATCH 2/6] update docs + bump version --- DESCRIPTION | 4 ++-- R/get_map_estimates.R | 7 ++++++- man/calc_residuals.Rd | 1 + man/get_map_estimates.Rd | 8 +++++++- man/ll_func_PKPDsim.Rd | 1 + man/ll_func_generic.Rd | 1 + 6 files changed, 18 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index d1f1985..9af98ba 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: PKPDmap Type: Package Title: MAP Bayesian estimates -Version: 1.0.3 -Date: 2025-06-18 +Version: 1.1.0 +Date: 2025-07-08 Author: Ron Keizer, Jasmine Hughes, Kara Woo Maintainer: Ron Keizer Description: Obtain MAP Bayesian estimates based on individual data and diff --git a/R/get_map_estimates.R b/R/get_map_estimates.R index 3c44605..d16d90c 100755 --- a/R/get_map_estimates.R +++ b/R/get_map_estimates.R @@ -39,7 +39,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` diff --git a/man/calc_residuals.Rd b/man/calc_residuals.Rd index 002eb5a..59b03e2 100644 --- a/man/calc_residuals.Rd +++ b/man/calc_residuals.Rd @@ -12,6 +12,7 @@ calc_residuals( parameters_population, covariates, regimen, + lagtime, omega_full, error, weights, diff --git a/man/get_map_estimates.Rd b/man/get_map_estimates.Rd index 02918ab..34b0cee 100755 --- a/man/get_map_estimates.Rd +++ b/man/get_map_estimates.Rd @@ -23,6 +23,7 @@ get_map_estimates( include_omega = TRUE, include_error = TRUE, regimen = NULL, + lagtime = NULL, t_init = 0, int_step_size = 0.1, ll_func = ll_func_PKPDsim, @@ -102,7 +103,12 @@ in this package.} \item{optimizer}{optimization library to use, default is `optim`} -\item{method}{optimization method, default `BFGS`} +\item{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.} \item{control}{list of options passed to `optim()` function} diff --git a/man/ll_func_PKPDsim.Rd b/man/ll_func_PKPDsim.Rd index b3aa654..74c3e5e 100644 --- a/man/ll_func_PKPDsim.Rd +++ b/man/ll_func_PKPDsim.Rd @@ -49,6 +49,7 @@ ll_func_PKPDsim( t_init = 0, calc_ofv, regimen, + lagtime = c(0), steady_state_analytic, include_omega, include_error, diff --git a/man/ll_func_generic.Rd b/man/ll_func_generic.Rd index 5d34f78..f5f2554 100644 --- a/man/ll_func_generic.Rd +++ b/man/ll_func_generic.Rd @@ -33,6 +33,7 @@ ll_func_generic( parameters, covariates = NULL, regimen = regimen, + lagtime = NULL, omega_full = omega_full, error = error, model, From 43826f491cedf18ce35dd0f3c6efaf8a8a495c21 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 11 Jul 2025 10:32:29 +0000 Subject: [PATCH 3/6] safer parsing obs_type --- DESCRIPTION | 4 ++-- R/parse_input_data.R | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9af98ba..797ae46 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: PKPDmap Type: Package Title: MAP Bayesian estimates -Version: 1.1.0 -Date: 2025-07-08 +Version: 1.1.1 +Date: 2025-07-11 Author: Ron Keizer, Jasmine Hughes, Kara Woo Maintainer: Ron Keizer Description: Obtain MAP Bayesian estimates based on individual data and diff --git a/R/parse_input_data.R b/R/parse_input_data.R index 62bd2fd..2267264 100644 --- a/R/parse_input_data.R +++ b/R/parse_input_data.R @@ -19,7 +19,7 @@ parse_input_data <- function( } } ## Parse data to include label for obs_type - if(is.null(obs_type_label)) { + if(is.null(obs_type_label) || is.null(data[[obs_type_label]])) { data$obs_type <- 1 } else { data$obs_type <- data[[obs_type_label]] From 97054693063fb8cba37b721795ded24c46f19e82 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 11 Jul 2025 10:49:22 +0000 Subject: [PATCH 4/6] fix missing parameter --- DESCRIPTION | 2 +- R/get_map_estimates.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 797ae46..bff5a46 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: PKPDmap Type: Package Title: MAP Bayesian estimates -Version: 1.1.1 +Version: 1.1.2 Date: 2025-07-11 Author: Ron Keizer, Jasmine Hughes, Kara Woo Maintainer: Ron Keizer diff --git a/R/get_map_estimates.R b/R/get_map_estimates.R index d16d90c..0862cee 100755 --- a/R/get_map_estimates.R +++ b/R/get_map_estimates.R @@ -145,7 +145,7 @@ get_map_estimates <- function( } ## Make sure lagtime is properly formatted (potentially pulled from model) - lagtime <- PKPDsim:::parse_lagtime(lagtime, model) + lagtime <- PKPDsim:::parse_lagtime(lagtime, model, parameters) ## Misc checks: check_inputs(model, data, parameters, omega, regimen, censoring, type) From 85f061dc88846461e0e4971e656fe4bd9e208015 Mon Sep 17 00:00:00 2001 From: Ron Keizer Date: Fri, 11 Jul 2025 10:59:41 +0000 Subject: [PATCH 5/6] clean up docs --- DESCRIPTION | 2 +- R/get_individual_parameters_from_fit.R | 3 ++- R/get_map_estimates.R | 2 ++ R/ll_func_PKPDsim.R | 12 +-------- R/ll_func_generic.R | 28 +++++++++++++++++++++ R/parse_weights.R | 1 + man/calc_residuals.Rd | 3 +++ man/get_individual_parameters_from_fit.Rd | 6 ++++- man/get_map_estimates.Rd | 3 +++ man/ll_func_PKPDsim.Rd | 30 ++++++++++++++--------- man/ll_func_generic.Rd | 18 ++++++++------ man/parse_weights.Rd | 2 ++ 12 files changed, 78 insertions(+), 32 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index bff5a46..69336a7 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: PKPDmap Type: Package Title: MAP Bayesian estimates -Version: 1.1.2 +Version: 1.1.3 Date: 2025-07-11 Author: Ron Keizer, Jasmine Hughes, Kara Woo Maintainer: Ron Keizer diff --git a/R/get_individual_parameters_from_fit.R b/R/get_individual_parameters_from_fit.R index 25c7a8e..8be8792 100644 --- a/R/get_individual_parameters_from_fit.R +++ b/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, diff --git a/R/get_map_estimates.R b/R/get_map_estimates.R index 0862cee..696c8ec 100755 --- a/R/get_map_estimates.R +++ b/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. diff --git a/R/ll_func_PKPDsim.R b/R/ll_func_PKPDsim.R index e595bf0..cd7f26e 100644 --- a/R/ll_func_PKPDsim.R +++ b/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 Date: Sat, 12 Jul 2025 07:06:58 +0000 Subject: [PATCH 6/6] fix parsing of data --- DESCRIPTION | 4 ++-- R/parse_input_data.R | 12 ++++++------ tests/testthat/test-parse_input_data.R | 5 +++-- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 69336a7..7e07728 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: PKPDmap Type: Package Title: MAP Bayesian estimates -Version: 1.1.3 -Date: 2025-07-11 +Version: 1.1.4 +Date: 2025-07-12 Author: Ron Keizer, Jasmine Hughes, Kara Woo Maintainer: Ron Keizer Description: Obtain MAP Bayesian estimates based on individual data and diff --git a/R/parse_input_data.R b/R/parse_input_data.R index 2267264..aefe23b 100644 --- a/R/parse_input_data.R +++ b/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) || is.null(data[[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,] } diff --git a/tests/testthat/test-parse_input_data.R b/tests/testthat/test-parse_input_data.R index 2457a2f..05b279e 100644 --- a/tests/testthat/test-parse_input_data.R +++ b/tests/testthat/test-parse_input_data.R @@ -9,11 +9,11 @@ test_that("parse_input_data handles regular data frame input", { result <- parse_input_data(data) # Check column names are lowercase - expect_equal(names(result), c("t", "obs_type", "y")) + expect_equal(names(result), c("t", "obs_type", "y", "obs_type")) # Check sorting expect_equal(result$t, c(1, 2, 3)) - expect_equal(result$obs_type, c(1, 1, 1)) + expect_equal(result$obs_type, c(1, 2, 1)) }) test_that("parse_input_data handles PKPDsim object", { @@ -133,3 +133,4 @@ test_that("parse_input_data preserves data values after transformations", { expect_equal(result$dv, c(10, 30)) expect_equal(result$t, c(2, 3)) }) +