diff --git a/NEWS.md b/NEWS.md index 1c8348f5..62abd419 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,6 +6,19 @@ the dosing including dose amount and route. # Development version +* Bug fix: `superposition()` no longer loops forever when steady-state is + requested (the default `n.tau=Inf`) and concentrations are extrapolated as + zero after `tlast` (for example, `auc.type="AUClast"` when `tlast < tau`). + Steady-state is now assessed on the concentrations that can become nonzero, + a warning is given when zero concentrations remain in the steady-state + profile, and an informative error (instead of an infinite loop) is raised if + steady-state cannot be reached within 10000 dosing intervals (#580). + Warnings, messages, and printed output from `superposition()` on a + `PKNCAconc` object are now raised by the parent process. Each subject runs + through `parallel::mclapply()`, which forks on every platform except + Windows, and a condition raised in a forked worker never reaches the + caller, so the new warning was previously invisible except on Windows. + * Corrected registry description strings that did not match the implementation (#582). `adj.r.squared.factor` is described as selecting the regression with the most points among those within the tolerance of diff --git a/R/superposition.R b/R/superposition.R index 32a1b96e..134e3468 100644 --- a/R/superposition.R +++ b/R/superposition.R @@ -35,6 +35,12 @@ #' times: 0 (zero), `dose.times`, `time modulo tau` (shifting `time` for each #' dose time as well), `additional.times`, and `tau`. #' +#' When concentrations are extrapolated as zero after `tlast` (for example, +#' `auc.type = "AUClast"` when `tlast < tau`), some concentrations in the +#' profile can never become nonzero. Steady-state (`n.tau = Inf`) is then +#' assessed on the nonzero concentrations, and a warning is given that zero +#' concentrations remain in the steady-state profile. +#' #' @seealso [interp.extrap.conc()] #' @family Superposition #' @export @@ -48,19 +54,38 @@ superposition.PKNCAconc <- function(conc, ...) { # Split the data by grouping and extract just the concentration and # time columns nested_data <- prepare_PKNCAconc(conc) + # parallel::mclapply() forks on every platform except Windows, and a + # condition raised in a forked worker never reaches the parent, so collect + # each worker's conditions with the result and re-emit them below. + # purrr::quietly() also captures messages and printed output, so those are + # re-emitted as well; otherwise adding a message() to the calculation later + # would silently discard it on every platform. + quiet_superposition <- purrr::quietly(superposition.numeric) tmp_results <- parallel::mclapply( X=seq_len(nrow(nested_data)), FUN=function(idx) { - superposition.numeric( + quiet_superposition( conc=nested_data$data_conc[[idx]]$conc, time=nested_data$data_conc[[idx]]$time, ... ) } ) + for (current_warning in unlist(lapply(tmp_results, FUN="[[", "warnings"))) { + warning(current_warning, call.=FALSE) + } + # `messages` keep their trailing newline; `output` has it stripped. + for (current_message in unlist(lapply(tmp_results, FUN="[[", "messages"))) { + message(current_message, appendLF=FALSE) + } + current_output <- unlist(lapply(tmp_results, FUN="[[", "output")) + current_output <- current_output[nzchar(current_output)] + if (length(current_output) > 0) { + writeLines(current_output) + } # Replace the concentration data with the new results - nested_data$data_conc <- tmp_results + nested_data$data_conc <- lapply(tmp_results, FUN="[[", "result") tidyr::unnest(nested_data, cols="data_conc") } @@ -212,10 +237,22 @@ superposition.numeric <- function(conc, time, dose.input = NULL, # Do the math! (Finally) current.tol <- steady.state.tol + 1 tau.count <- 0 + # Generous backstop so that a case that cannot reach steady-state errors + # instead of looping forever (normal convergence takes far fewer + # intervals; the default steady.state.tol needs at most ~1000 even for + # negligible accumulation decay). + max.tau.count <- 10000L # Stop either for reaching steady-state or for reaching the requested number of doses while (tau.count < n.tau & !is.na(current.tol) & current.tol >= steady.state.tol) { + if (is.infinite(n.tau) && tau.count >= max.tau.count) { + stop( + "Superposition did not reach steady-state within ", max.tau.count, + " dosing intervals (steady.state.tol=", steady.state.tol, + "). Check the lambda.z, tau, and steady.state.tol inputs, or use a finite n.tau." + ) + } prev.conc <- ret$conc # Perform the dosing for a single dosing interval. for (i in seq_along(dose.times)) { @@ -238,14 +275,29 @@ superposition.numeric <- function(conc, time, dose.input = NULL, ) } tau.count <- tau.count + 1 - if (any(ret$conc %in% 0)) { - # prevent division by 0. Since not all concentrations are 0, - # all values will eventually be nonzero. - current.tol <- steady.state.tol + 1 - } else { + zero.conc <- ret$conc %in% 0 + if (!any(zero.conc)) { current.tol <- max(1-(prev.conc/ret$conc)) + } else if ((tau*(tau.count - 1) - max(dose.times)) > tlast) { + # The dosing interval just added drew all of its concentrations from + # times after tlast, where extrapolation may be zero (e.g. with + # auc.type="AUClast"). A concentration that is still zero can + # therefore never become nonzero (a structural zero), so assess + # steady-state on the nonzero concentrations only. + current.tol <- max(1-(prev.conc[!zero.conc]/ret$conc[!zero.conc])) + } else { + # A zero concentration may still become nonzero from a later dosing + # interval; prevent division by 0 and do not yet assess steady-state. + current.tol <- steady.state.tol + 1 } } + if (!is.na(current.tol) && current.tol < steady.state.tol && any(ret$conc %in% 0)) { + warning( + "Zero concentrations remain in the steady-state superposition profile. ", + "They come from concentrations extrapolated as zero after tlast (", + signif(tlast, 6), ") with auc.type=", dQuote(auc.type, q = FALSE), "." + ) + } } ret } diff --git a/man/superposition.Rd b/man/superposition.Rd index cf427bcc..6c59aa28 100644 --- a/man/superposition.Rd +++ b/man/superposition.Rd @@ -93,6 +93,12 @@ Compute noncompartmental superposition for repeated dosing The returned superposition times will include all of the following times: 0 (zero), \code{dose.times}, \verb{time modulo tau} (shifting \code{time} for each dose time as well), \code{additional.times}, and \code{tau}. + +When concentrations are extrapolated as zero after \code{tlast} (for example, +\code{auc.type = "AUClast"} when \code{tlast < tau}), some concentrations in the +profile can never become nonzero. Steady-state (\code{n.tau = Inf}) is then +assessed on the nonzero concentrations, and a warning is given that zero +concentrations remain in the steady-state profile. } \seealso{ \code{\link[=interp.extrap.conc]{interp.extrap.conc()}} diff --git a/tests/testthat/test-superpostion.R b/tests/testthat/test-superpostion.R index 15c9c340..e9313605 100644 --- a/tests/testthat/test-superpostion.R +++ b/tests/testthat/test-superpostion.R @@ -514,3 +514,154 @@ test_that("PKNCAconc superposition", { ) ) }) + +test_that("superposition reaches steady-state despite structural zero concentrations (issue 580)", { + # Theoph subjects 6 and 10 have tlast < tau=24, so with auc.type="AUClast" + # the concentrations after tlast extrapolate as exactly zero and stay zero at + # steady-state. This formerly looped forever with n.tau=Inf because the + # zero-concentration guard reset the steady-state tolerance every interval. + d_theoph_6 <- datasets::Theoph[datasets::Theoph$Subject == 6, ] + expect_warning( + v_ss <- + superposition( + conc=d_theoph_6$conc, time=d_theoph_6$Time, tau=24, auc.type="AUClast" + ), + regexp='Zero concentrations remain in the steady-state superposition profile. They come from concentrations extrapolated as zero after tlast (23.85) with auc.type="AUClast".', + fixed=TRUE + ) + # No accumulation is possible because every concentration one or more + # dosing intervals after dosing is zero, so the steady-state profile is + # exactly the single-dose profile with structural zeros at times 0 and tau. + expect_identical( + v_ss, + data.frame(conc=c(d_theoph_6$conc, 0), time=c(d_theoph_6$Time, 24)) + ) + # A finite n.tau with the same input gave the same concentrations before the + # fix (all added dosing intervals contribute zero); that is preserved, and + # partial accumulation with a finite n.tau does not warn. + expect_silent( + v_1 <- + superposition( + conc=d_theoph_6$conc, time=d_theoph_6$Time, tau=24, n.tau=1, auc.type="AUClast" + ) + ) + expect_identical(v_ss, v_1) + + # The default auc.type="AUCinf" on the same subject is unaffected by the fix: + # all concentrations become positive and steady-state is reached silently + # (values are pinned from before the fix). + expect_silent( + v_aucinf <- superposition(conc=d_theoph_6$conc, time=d_theoph_6$Time, tau=24) + ) + expect_equal( + v_aucinf, + data.frame( + conc=c(1.03361737415698, 2.29940375348928, 4.06230162309106, + 7.37435349523853, 7.18488330467035, 6.28550707796221, + 5.60636744789876, 4.57905606419521, 3.92005369868033, + 3.13726988893256, 1.04734393070531, 1.03364150451656), + time=c(0, 0.27, 0.58, 1.15, 2.03, 3.57, 5, 7, 9.22, 12.1, 23.85, 24) + ) + ) + + # AUClast with tlast > tau converges without a warning: the time-zero + # concentration is zero after the first dosing interval, becomes nonzero on + # the second, and no zeros remain at steady-state. + d_theoph_1 <- datasets::Theoph[datasets::Theoph$Subject == 1, ] + expect_silent( + v_s1 <- + superposition( + conc=d_theoph_1$conc, time=d_theoph_1$Time, tau=24, auc.type="AUClast", + check.blq=FALSE + ) + ) + expect_true(all(v_s1$conc > 0)) + + # The PKNCAconc method (the form reported in issue 580) returns with one + # warning per affected subject (check.blq=FALSE because subject 10 has a + # nonzero concentration at time 0) + d_theoph_6_10 <- + datasets::Theoph[datasets::Theoph$Subject %in% c(6, 10), ] + conc_obj <- PKNCAconc(conc~Time|Subject, data=d_theoph_6_10) + w_obj <- + capture_warnings( + v_obj <- superposition(conc_obj, tau=24, auc.type="AUClast", check.blq=FALSE) + ) + expect_length(w_obj, 2) + expect_match( + w_obj, + regexp="Zero concentrations remain in the steady-state superposition profile", + fixed=TRUE, + all=TRUE + ) + expect_true(all(is.finite(v_obj$conc))) +}) + +test_that("superposition errors instead of hanging when steady-state cannot be reached", { + # A negligible lambda.z makes each added dosing interval contribute ~clast, + # so the relative change shrinks only like 1/n.tau and steady.state.tol=1e-8 + # cannot be reached within the iteration cap. + expect_error( + superposition( + conc=c(0, 1, 0.5), time=0:2, tau=1, + lambda.z=1e-8, clast.pred=0.5, tlast=2, + steady.state.tol=1e-8 + ), + regexp="Superposition did not reach steady-state within 10000 dosing intervals", + fixed=TRUE + ) +}) + +test_that("superposition warns about structural zeros for a single profile (issue 580)", { + # superposition.PKNCAconc() runs each subject through parallel::mclapply(), + # which forks everywhere except Windows, and a warning raised in a forked + # worker never reaches the parent; the method collects them with + # purrr::quietly() and re-emits them. The PKNCAconc test above cannot show + # that on Windows, where mclapply() falls back to lapply(), so pin the + # underlying warning here where it is raised. + d_theoph_6 <- datasets::Theoph[datasets::Theoph$Subject == 6, ] + expect_warning( + v_6 <- + superposition( + conc=d_theoph_6$conc, time=d_theoph_6$Time, + tau=24, auc.type="AUClast", check.blq=FALSE + ), + regexp="Zero concentrations remain in the steady-state superposition profile", + fixed=TRUE + ) + expect_equal(nrow(v_6), 12) + expect_true(all(is.finite(v_6$conc))) +}) + +test_that("superposition.PKNCAconc re-emits everything its workers produce (issue 580)", { + # parallel::mclapply() forks everywhere except Windows, and conditions raised + # in a forked worker never reach the parent. purrr::quietly() collects the + # warnings, messages, and printed output so the method can re-emit them. + # Mock the calculation so all three are produced regardless of the data. + d_theoph_6_10 <- + datasets::Theoph[datasets::Theoph$Subject %in% c(6, 10), ] + conc_obj <- PKNCAconc(conc~Time|Subject, data=d_theoph_6_10) + local_mocked_bindings( + superposition.numeric=function(conc, time, ...) { + warning("mocked warning") + message("mocked message") + cat("mocked output", fill=TRUE) + data.frame(conc=c(1, 2), time=c(0, 1)) + } + ) + printed <- NULL + warned <- + capture_warnings( + messaged <- + capture_messages( + # Assign the result so capture.output() does not also see it printed. + printed <- capture.output(v_obj <- superposition(conc_obj, tau=24)) + ) + ) + # One of each per subject, and messages are not given an extra blank line. + expect_length(warned, 2) + expect_equal(unique(warned), "mocked warning") + expect_length(messaged, 2) + expect_equal(unique(trimws(messaged)), "mocked message") + expect_equal(unique(printed), "mocked output") +})