Add exclude_nca (Fix #20)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: PKNCA
 Type: Package
 Title: Perform Pharmacokinetic Non-Compartmental Analysis
-Version: 0.8.1.9016
+Version: 0.8.1.9017
 Date: 2017-09-09
 Authors@R: c(
     person("Bill", "Denney", email="wdenney@humanpredictions.com", role=c("aut", "cre")),

NAMESPACE

@@ -74,6 +74,9 @@ export(choose.auc.intervals)
 export(clean.conc.blq)
 export(clean.conc.na)
 export(exclude)
+export(exclude_nca_max.aucinf.pext)
+export(exclude_nca_min.hl.r.squared)
+export(exclude_nca_span.ratio)
 export(extrapolate.conc)
 export(find.tau)
 export(findOperator)
@@ -158,6 +161,7 @@ importFrom(nlme,nlme)
 importFrom(nlme,ranef)
 importFrom(parallel,mclapply)
 importFrom(plyr,rbind.fill)
+importFrom(rlang,syms)
 importFrom(stats,AIC)
 importFrom(stats,as.formula)
 importFrom(stats,coef)

NEWS.md

@@ -20,6 +20,7 @@ the dosing including dose amount and route.
 * BACKWARD INCOMPATIBILITY: Functions that previously may have
   returned Infinity due to dividing by zero (e.g. when AUC=0
   calculating clearance) now return NA.
+* exclude_nca* functions added (Fixes issue #20)
 * Add manual half-life point selection (Fixes issue #18)
 * Add parameters for Ceoi and intravenous MRT
 * Updated vignettes to improve clarity

R/exclude.R

@@ -31,6 +31,7 @@
 #'         mask=c(TRUE, rep(FALSE, 6)))
 #' @export
 #' @importFrom dplyr "%>%"
+#' @importFrom rlang syms
 exclude <- function(object, reason, mask, FUN)
   UseMethod("exclude")
 

R/exclude_nca.R

@@ -0,0 +1,116 @@
+#' Exclude NCA parameters based on examining the parameter set.
+#' 
+#' @param min.span.ratio The minimum acceptable span ratio (uses 
+#'   \code{PKNCA.options("min.span.ratio")} if not provided).
+#' @param max.aucinf.pext The maximum acceptable percent AUC
+#'   extrapolation (uses \code{PKNCA.options("max.aucinf.pext")} if not
+#'   provided).
+#' @param min.hl.r.squared The minimum acceptable r-squared value for
+#'   half-life (uses \code{PKNCA.options("min.hl.r.squared")} if not
+#'   provided).
+#' @examples
+#' my_conc <- PKNCAconc(data.frame(conc=1.1^(3:0), time=0:3, subject=1), conc~time|subject)
+#' my_data <- PKNCAdata(my_conc, intervals=data.frame(start=0, end=Inf, aucinf.obs=TRUE, aucpext.obs=TRUE))
+#' my_result <- pk.nca(my_data)
+#' my_result_excluded <- exclude(my_result, FUN=exclude_nca_max.aucinf.pext())
+#' as.data.frame(my_result_excluded)
+#' @name exclude_nca
+NULL
+
+#' @describeIn exclude_nca Exclude based on span.ratio
+#' @export
+exclude_nca_span.ratio <- function(min.span.ratio) {
+  affected_parameters <- get.parameter.deps("half.life")
+  missing_min.span.ratio <- missing(min.span.ratio)
+  function(x, ...) {
+    if (missing_min.span.ratio) {
+      min.span.ratio <- PKNCA.options("min.span.ratio")
+    }
+    ret <- rep(NA_character_, nrow(x))
+    if (!is.na(min.span.ratio)) {
+      idx_span.ratio <- which(x$PPTESTCD %in% "span.ratio")
+      if (length(idx_span.ratio) == 0) {
+        # Do nothing, it wasn't calculated
+      } else if (length(idx_span.ratio) == 1) {
+        current_span.ratio <- x$PPORRES[idx_span.ratio]
+        drop_span_ratio <-
+          !is.na(current_span.ratio) &
+          current_span.ratio < min.span.ratio
+        if (drop_span_ratio) {
+          ret[x$PPTESTCD %in% affected_parameters] <-
+            sprintf("Span ratio < %g", min.span.ratio)
+        }
+      } else if (length(idx_span.ratio) > 1) {
+        stop("Should not see more than one span.ratio")
+      }
+    }
+    ret
+  }
+}
+
+#' @describeIn exclude_nca Exclude based on AUC percent extrapolated
+#'   (both observed and predicted)
+#' @export
+exclude_nca_max.aucinf.pext <-  function(max.aucinf.pext) {
+  affected_parameters <-
+    list(obs=get.parameter.deps("aucinf.obs"),
+         pred=get.parameter.deps("aucinf.pred"))
+  missing_max.aucinf.pext <- missing(max.aucinf.pext)
+  function(x, ...) {
+    if (missing_max.aucinf.pext) {
+      max.aucinf.pext <- PKNCA.options("max.aucinf.pext")
+    }
+    ret <- rep(NA_character_, nrow(x))
+    if (!is.na(max.aucinf.pext)) {
+      for (ext_type in c("obs", "pred")) {
+        idx_pext <- which(x$PPTESTCD %in% paste0("aucpext.", ext_type))
+        if (length(idx_pext) == 0) {
+          # Do nothing, it wasn't calculated
+        } else if (length(idx_pext) == 1) {
+          current_pext <- x$PPORRES[idx_pext]
+          drop_pext <-
+            !is.na(current_pext) &
+            current_pext > max.aucinf.pext
+          if (drop_pext) {
+            ret[x$PPTESTCD %in% affected_parameters[[ext_type]]] <-
+              sprintf("AUC percent extrapolated > %g", max.aucinf.pext)
+          }
+        } else if (length(idx_pext) > 1) {
+          stop("Should not see more than one aucpext.", ext_type)
+        }
+      }
+    }
+    ret
+  }
+}
+
+#' @describeIn exclude_nca Exclude based on half-life r-squared
+#' @export
+exclude_nca_min.hl.r.squared <- function(min.hl.r.squared) {
+  affected_parameters <- get.parameter.deps("half.life")
+  missing_min.hl.r.squared <- missing(min.hl.r.squared)
+  function(x, ...) {
+    if (missing_min.hl.r.squared) {
+      min.hl.r.squared <- PKNCA.options("min.hl.r.squared")
+    }
+    ret <- rep(NA_character_, nrow(x))
+    if (!is.na(min.hl.r.squared)) {
+      idx_r.squared <- which(x$PPTESTCD %in% "r.squared")
+      if (length(idx_r.squared) == 0) {
+        # Do nothing, it wasn't calculated
+      } else if (length(idx_r.squared) == 1) {
+        current_r.squared <- x$PPORRES[idx_r.squared]
+        drop_r.squared <-
+          !is.na(current_r.squared) &
+          current_r.squared < min.hl.r.squared
+        if (drop_r.squared) {
+          ret[x$PPTESTCD %in% affected_parameters] <-
+            sprintf("Half-life r-squared < %g", min.hl.r.squared)
+        }
+      } else if (length(idx_r.squared) > 1) {
+        stop("Should not see more than one r.squared")
+      }
+    }
+    ret
+  }
+}

tests/testthat/test-exclude_nca.R

@@ -0,0 +1,47 @@
+context("exclude_nca")
+
+test_that("exclude_nca", {
+  my_conc <- PKNCAconc(data.frame(conc=c(1.1^(3:0), 1.1), time=0:4, subject=1), conc~time|subject)
+  my_data <- PKNCAdata(my_conc, intervals=data.frame(start=0, end=Inf, aucinf.obs=TRUE, aucpext.obs=TRUE))
+  my_result <- pk.nca(my_data)
+
+  my_result_excluded <- exclude(my_result, FUN=exclude_nca_max.aucinf.pext())
+  expect_equal(as.data.frame(my_result_excluded)$exclude,
+               c(rep(NA_character_, nrow(my_result_excluded$result)-2),
+                 rep("AUC percent extrapolated > 20", 2)))
+
+  my_result_excluded <- exclude(my_result, FUN=exclude_nca_max.aucinf.pext(50))
+  expect_equal(as.data.frame(my_result_excluded)$exclude,
+               c(rep(NA_character_, nrow(my_result_excluded$result)-2),
+                 rep("AUC percent extrapolated > 50", 2)))
+
+  my_result_excluded <- exclude(my_result, FUN=exclude_nca_span.ratio())
+  expect_equal(as.data.frame(my_result_excluded)$exclude,
+               c(rep(NA_character_, 4),
+                 rep("Span ratio < 2", 10)))
+  my_result_excluded <- exclude(my_result, FUN=exclude_nca_span.ratio(1))
+  expect_equal(as.data.frame(my_result_excluded)$exclude,
+               c(rep(NA_character_, 4),
+                 rep("Span ratio < 1", 10)))
+
+  my_result_excluded <- exclude(my_result, FUN=exclude_nca_min.hl.r.squared())
+  expect_equal(as.data.frame(my_result_excluded)$exclude,
+               c(rep(NA_character_, 4),
+                 rep("Half-life r-squared < 0.9", 10)))
+  my_result_excluded <- exclude(my_result, FUN=exclude_nca_min.hl.r.squared(0.95))
+  expect_equal(as.data.frame(my_result_excluded)$exclude,
+               c(rep(NA_character_, 4),
+                 rep("Half-life r-squared < 0.95", 10)))
+
+  my_data <- PKNCAdata(my_conc, intervals=data.frame(start=0, end=Inf, cmax=TRUE))
+  my_result <- pk.nca(my_data)
+  expect_equal(my_result,
+               exclude(my_result, FUN=exclude_nca_max.aucinf.pext()),
+               info="Result is ignored when not calculated")
+  expect_equal(my_result,
+               exclude(my_result, FUN=exclude_nca_span.ratio()),
+               info="Result is ignored when not calculated")
+  expect_equal(my_result,
+               exclude(my_result, FUN=exclude_nca_min.hl.r.squared()),
+               info="Result is ignored when not calculated")
+})