Merge pull request #13 from nlmixr2/12-allow-vpcplot-to-use-vpc-simul…

…ation-from-vpcsim Take care of vpcSim objects
nlmixr2 · Sep 6, 2022 · 55986f0 · 55986f0
2 parents 91a6cab + da6452b
commit 55986f0
Show file tree

Hide file tree

Showing 5 changed files with 131 additions and 31 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -2,12 +2,12 @@ Package: nlmixr2plot
 Title: Nonlinear Mixed Effects Models in Population PK/PD, Plot Functions
 Version: 2.0.6
 Authors@R: c(
-    person("Matthew", "Fidler", , "matthew.fidler@gmail.com", role = c("aut", "cre"),
+    person("Matthew", "Fidler", email="matthew.fidler@gmail.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-8538-6691")),
-    person("Bill", "Denney", , "wdenney@humanpredictions.com", role = "ctb",
+    person("Bill", "Denney", email="wdenney@humanpredictions.com", role = "ctb",
            comment = c(ORCID = "0000-0002-5759-428X")),
-    person("Wenping", "Wang", , "wwang8198@gmail.com", role = "aut"),
-    person("Vipul", "Mann", , "vm2583@columbia.edu", role = "aut")
+    person("Wenping", "Wang", email="wwang8198@gmail.com", role = "aut"),
+    person("Vipul", "Mann", email="vm2583@columbia.edu", role = "aut")
   )
 Description: Fit and compare nonlinear mixed-effects models in
     differential equations with flexible dosing information commonly seen
@@ -35,4 +35,4 @@ Suggests:
 Config/testthat/edition: 3
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.1
diff --git a/NAMESPACE b/NAMESPACE
@@ -13,7 +13,10 @@ export(geom_cens)
 export(stat_amt)
 export(stat_cens)
 export(traceplot)
+export(vpcCens)
+export(vpcCensTad)
 export(vpcPlot)
+export(vpcPlotTad)
 importFrom(ggplot2,.data)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,element_blank)

diff --git a/R/vpcPlot.R b/R/vpcPlot.R
@@ -2,10 +2,17 @@
 #'
 #' @param fit nlmixr2 fit object
 #' @param data this is the data to use to augment the VPC fit.  By
-#'     default is the fitted data, (can be retrieved by
-#'     \code{\link[nlme]{getData}}), but it can be changed by specifying
-#'     this argument.
+#'   default is the fitted data, (can be retrieved by
+#'   \code{\link[nlme]{getData}}), but it can be changed by specifying
+#'   this argument.
 #' @param n Number of VPC simulations.  By default 100
+#' @param idv Name of independent variable. For `vpcPlot()` and
+#'   `vpcCens()` the default is `"time"` for `vpcPlotTad()` and
+#'   `vpcCensTad()` this is `"tad"`
+#' @param cens is a boolean to show if this is a censoring plot or
+#'   not.  When `cens=TRUE` this is actually a censoring vpc plot
+#'   (with `vpcCens()` and `vpcCensTad()`).  When `cens=FALSE` this is
+#'   traditional VPC plot (`vpcPlot()` and `vpcPlotTad()`).
 #' @inheritParams vpc::vpc
 #' @inheritParams rxode2::rxSolve
 #' @param ... Args sent to \code{\link[rxode2]{rxSolve}}
@@ -48,13 +55,16 @@ vpcPlot <- function(fit, data = NULL, n = 300, bins = "jenks",
                     n_bins = "auto", bin_mid = "mean",
                     show = NULL, stratify = NULL, pred_corr = FALSE,
                     pred_corr_lower_bnd = 0, pi = c(0.05, 0.95), ci = c(0.05, 0.95),
-                    uloq = NULL, lloq = NULL, log_y = FALSE, log_y_min = 0.001,
+                    uloq = fit$dataUloq, lloq = fit$dataLloq, log_y = FALSE, log_y_min = 0.001,
                     xlab = NULL, ylab = NULL, title = NULL, smooth = TRUE, vpc_theme = NULL,
                     facet = "wrap", scales = "fixed", labeller = NULL, vpcdb = FALSE,
-                    verbose = FALSE, ..., seed=1009) {
+                    verbose = FALSE, ..., seed=1009,
+                    idv="time", cens=FALSE) {
+  force(idv)
   rxode2::rxReq("vpc")
   .ui <- fit$ui
-  .obsLst <- .vpcUiSetupObservationData(fit, data)
+  .obsLst <- .vpcUiSetupObservationData(fit, data=data, idv=idv)
+  .obs <- .obsLst$obs
   .no <- .obsLst$namesObs
   .nol <- .obsLst$namesObsLower
   .obs <- .obsLst$obs
@@ -72,12 +82,16 @@ vpcPlot <- function(fit, data = NULL, n = 300, bins = "jenks",
     }
   }
   # Simulate with VPC
-  .sim <- nlmixr2est::vpcSim(fit, ..., keep=stratify, n=n, pred=pred_corr, seed=seed)
-  .sim <- nlmixr2est::vpcSimExpand(fit, .sim, stratify)
+  if (inherits(fit, "nlmixr2vpcSim")) {
+    .sim <- fit
+  } else {
+    .sim <- nlmixr2est::vpcSim(fit, ..., keep=stratify, n=n, pred=pred_corr, seed=seed)
+  }
+  .sim <- nlmixr2est::vpcSimExpand(fit, .sim, stratify, .obs)
   .simCols <- list(
     id="id",
     dv="sim",
-    idv="time")
+    idv=idv)
   if (pred_corr) {
     .simCols <- c(.simCols, list(pred="pred"))
     .si <- nlmixr2est::.nlmixr2estLastPredSimulationInfo()
@@ -99,7 +113,7 @@ vpcPlot <- function(fit, data = NULL, n = 300, bins = "jenks",
     .w <- which(.no == "sim")
     names(.obs)[.w] <- "pred"
     .obsCols$pred <- "pred"
-    .obsCols$idv <- "time"
+    .obsCols$idv <- idv
     .obsCols$id <- "id"
     if (any(names(.obs) == "dv")) {
       .obsCols$dv <- "dv"
@@ -124,14 +138,50 @@ vpcPlot <- function(fit, data = NULL, n = 300, bins = "jenks",
   if (length(.w) > 0) {
     .obs <- .obs[, -.w]
   }
-  vpc::vpc_vpc(sim=.sim, sim_cols=.simCols,
-               obs=.obs, obs_cols=.obsCols,
-               bins=bins, n_bins=n_bins, bin_mid=bin_mid,
-               show = show, stratify = stratify, pred_corr = pred_corr,
-               pred_corr_lower_bnd = pred_corr_lower_bnd, pi = pi, ci = ci,
-               uloq = uloq, lloq = lloq, log_y = log_y, log_y_min = log_y_min,
-               xlab = xlab, ylab = ylab, title = title, smooth = smooth, vpc_theme = vpc_theme,
-               facet = facet, scales=scales, labeller = labeller, vpcdb = vpcdb, verbose = verbose)
+  .obsCols$idv <- idv
+  .w <- which(tolower(names(.sim)) == "id")
+  names(.sim)[.w] <- "id"
+  if (cens) {
+    if (is.null(lloq) && is.null(uloq)) {
+      stop("this data is not censored")
+    }
+    vpc::vpc_cens(sim=.sim, sim_cols=.simCols,
+                 obs=.obs, obs_cols=.obsCols,
+                 bins=bins, n_bins=n_bins, bin_mid=bin_mid,
+                 show = show, stratify = stratify, 
+                 ci = ci,
+                 uloq = uloq, lloq = lloq, 
+                 xlab = xlab, ylab = ylab, title = title, smooth = smooth, vpc_theme = vpc_theme,
+                 facet = facet, labeller = labeller, vpcdb = vpcdb, verbose = verbose)    
+  } else {
+    vpc::vpc_vpc(sim=.sim, sim_cols=.simCols,
+                 obs=.obs, obs_cols=.obsCols,
+                 bins=bins, n_bins=n_bins, bin_mid=bin_mid,
+                 show = show, stratify = stratify, pred_corr = pred_corr,
+                 pred_corr_lower_bnd = pred_corr_lower_bnd, pi = pi, ci = ci,
+                 uloq = uloq, lloq = lloq, log_y = log_y, log_y_min = log_y_min,
+                 xlab = xlab, ylab = ylab, title = title, smooth = smooth, vpc_theme = vpc_theme,
+                 facet = facet, scales=scales, labeller = labeller, vpcdb = vpcdb, verbose = verbose)
+  }
+}
+
+#' @rdname vpcPlot
+#' @export
+vpcPlotTad <- function(..., idv="tad") {
+  vpcPlot(..., idv=idv)
+}
+
+
+#' @rdname vpcPlot
+#' @export
+vpcCensTad <- function(..., cens=TRUE, idv="tad") {
+  vpcPlot(..., cens=cens, idv=idv)
+}
+
+#' @rdname vpcPlot
+#' @export
+vpcCens <- function(..., cens=TRUE, idv="time") {
+  vpcPlot(..., cens=cens, idv=idv)
 }
 
 #' Setup Observation data for VPC
@@ -141,7 +191,7 @@ vpcPlot <- function(fit, data = NULL, n = 300, bins = "jenks",
 #' @return List with `namesObs`, `namesObsLower`, `obs` and `obsCols`
 #' @author Matthew L. Fidler
 #' @noRd
-.vpcUiSetupObservationData <- function(fit, data=NULL) {
+.vpcUiSetupObservationData <- function(fit, data=NULL, idv="time") {
   if (!is.null(data)) {
     .obs <- data
   } else {
@@ -163,10 +213,23 @@ vpcPlot <- function(fit, data = NULL, n = 300, bins = "jenks",
   }
   .obsCols <- c(.obsCols,
                 list(dv=.no[.wo]))
-  .wo <- which(.nol == "time")
+  .wo <- which(.nol == idv)
   if (length(.wo) != 1) {
-    stop("cannot find 'time' in original dataset",
-         call.=FALSE)
+    if (any(names(fit) == idv)) {
+      .fit <- as.data.frame(fit)
+      .wid <- which(tolower(names(.fit)) == "id")
+      names(.fit)[.wid] <- "ID"
+      .fit$nlmixrRowNums <-  fit$env$.rownum
+      .fit <- .fit[, c("ID", idv, "nlmixrRowNums")]
+      .obs$nlmixrRowNums <- seq_along(.obs$ID)
+      .obs <- merge(.obs, .fit, by=c("ID", "nlmixrRowNums"), all.x=TRUE)
+      .wo <- which(.nol == idv)
+    } else {
+      stop("cannot find '", idv, "' in original dataset",
+           call.=FALSE)
+    }
+  } else {
+    names(.obs)[.wo] <- idv
   }
   .obsCols <- c(.obsCols,
                 list(idv=.no[.wo]))

diff --git a/man/vpcPlot.Rd b/man/vpcPlot.Rd
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
@@ -22,13 +22,27 @@ test_that("test plots with vdiffr", {
     })
   }
 
-  fit <- nlmixr2est::nlmixr(one.cmt, nlmixr2data::theo_sd, est="focei",
+  censData <- nlmixr2data::theo_md
+  # Assign CENS = 1 for bloq values, otherwise CENS = 0.
+  censData$CENS[censData$DV < 3 & censData$AMT == 0] <- 1
+  censData$CENS[censData$DV >= 3 & censData$AMT == 0] <- 0
+  #
+  # Set DV to LOQ for all censored items
+  censData$DV[censData$CENS == 1] <-  3
+#
+
+  fit <- nlmixr2est::nlmixr(one.cmt, censData,
+                            est="focei",
                             table=nlmixr2est::tableControl(npde=TRUE))
 
   apo <- nlmixr2est::augPred(fit)
   expect_error(plot(apo), NA)
   expect_error(vpcPlot(fit), NA)
+  expect_error(vpcPlotTad(fit), NA)
   expect_error(vpcPlot(fit, pred_corr=TRUE), NA)
+  expect_error(vpcPlotTad(fit, pred_corr=TRUE), NA)
+  expect_error(vpcCens(fit), NA)
+  expect_error(vpcCensTad(fit), NA)
 
   expect_error(plot(fit), NA)