Merge pull request #23 from nlmixr2/nameplot

Make returned object hierarchical to give names to all objects
nlmixr2 · Nov 14, 2022 · 764d033 · 764d033
2 parents eed98e3 + 456eb36
commit 764d033
Show file tree

Hide file tree

Showing 5 changed files with 168 additions and 175 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,9 @@
+# nlmixr2plot 2.0.8
+
+* `plot()` now returns a named list of lists so that users can more easily
+  choose which plots to include, if all plots are not desired.  Or the user
+  could use those names as the basis of figure captions (fix #8).
+
 # nlmixr2plot 2.0.7
 
 * Added fixes for upcoming 'ggplot2' 3.4 release

diff --git a/R/plot.nlmixr2.R b/R/plot.nlmixr2.R
@@ -126,25 +126,23 @@
 #' Plot some standard goodness of fit plots for the focei fitted object
 #'
 #' @param x a focei fit object
-#' @param ... additional arguments
-#' @return Nothing, called for its side effects
+#' @param ... additional arguments (currently ignored)
+#' @return An `nlmixr2PlotList` object (a list of ggplot2 objects with easier
+#'   plotting for all of them at the same time)
 #' @author Wenping Wang & Matthew Fidler
 #' @examples
 #' \donttest{
-#'
 #' library(nlmixr2est)
-#' ## The basic model consiss of an ini block that has initial estimates
 #' one.compartment <- function() {
 #'   ini({
-#'     tka <- 0.45 # Log Ka
-#'     tcl <- 1 # Log Cl
-#'     tv <- 3.45    # Log V
+#'     tka <- 0.45
+#'     tcl <- 1
+#'     tv <- 3.45
 #'     eta.ka ~ 0.6
 #'     eta.cl ~ 0.3
 #'     eta.v ~ 0.1
 #'     add.sd <- 0.7
 #'   })
-#'   # and a model block with the error sppecification and model specification
 #'   model({
 #'     ka <- exp(tka + eta.ka)
 #'     cl <- exp(tcl + eta.cl)
@@ -157,153 +155,139 @@
 #' }
 #'
 #' ## The fit is performed by the function nlmixr/nlmix2 specifying the model, data and estimate
-#' fit <- nlmixr2(one.compartment, theo_sd,  est="saem", saemControl(print=0))
+#' fit <- nlmixr2(one.compartment, theo_sd,  est="saem", saemControl(print=0, nBurn = 10, nEm = 20))
 #'
 #' # This shows many goodness of fit plots
 #' plot(fit)
-#'
 #' }
 #' @export
 plot.nlmixr2FitData <- function(x, ...) {
   .lst <- list()
   object <- x
   .tp <- traceplot(x)
   if (!is.null(.tp)) {
-    .lst[[length(.lst) + 1]] <- .tp
+    .lst[["traceplot"]] <- .tp
   }
   if (exists(".bootPlotData", object$env)) {
     .bp <- nlmixr2extra::bootplot(x)
-    .lst[[length(.lst) + 1]] <- .bp
+    .lst[["bootplot"]] <- .bp
   }
   .dat <- .setupPlotData(x)
-  .hasCwres <- any(names(.dat) == "CWRES")
-  .hasNpde <- any(names(.dat) == "NPD")
-  .hasPred <- any(names(.dat) == "PRED")
-  .hasIpred <- any(names(.dat) == "IPRED")
   for (.cmt in levels(.dat$CMT)) {
-    .dat0 <- .dat[.dat$CMT == .cmt,, drop = FALSE]
-    if (dim(.dat0)[1] > 0) {
-      if (.hasPred) {
-        .p1 <- .dvPlot(.dat0, c("PRED", "IPRED")) +
-          ggplot2::ggtitle(.cmt, "DV vs PRED/IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
+    .lst[[.cmt]] <- plotCmt(.dat, cmt = .cmt)
+  }
 
-        .p1 <- .dvPlot(.dat0, c("PRED", "IPRED"), TRUE) +
-          ggplot2::ggtitle(.cmt, "log-scale DV vs PRED/IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
-      } else if (.hasIpred) {
-        .p1 <- .dvPlot(.dat0, "IPRED") +
-          ggplot2::ggtitle(.cmt, "DV vs IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
+  class(.lst) <- "nlmixr2PlotList"
+  .lst
+}
 
-        .p1 <- .dvPlot(.dat0, "IPRED", TRUE) +
-          ggplot2::ggtitle(.cmt, "log-scale DV vs IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
-      }
+#' Plot data from one compartment
+#'
+#' @inheritParams plot.nlmixr2FitData
+#' @param cmt The value of the current compartment
+#' @return A list of ggplot2 objects
+#' @noRd
+plotCmt <- function(x, cmt) {
+  .lst <- list()
+  .hasCwres <- any(names(x) == "CWRES")
+  .hasNpde <- any(names(x) == "NPD")
+  .hasPred <- any(names(x) == "PRED")
+  .hasIpred <- any(names(x) == "IPRED")
+  .datCmt <- x[x$CMT == cmt,, drop = FALSE]
+  if (nrow(.datCmt) > 0) {
+    if (.hasPred) {
+      .lst[["dv_pred_ipred_linear"]] <-
+        .dvPlot(.datCmt, c("PRED", "IPRED")) +
+        ggplot2::ggtitle(cmt, "DV vs PRED/IPRED")
 
+      .lst[["dv_pred_ipred_log"]] <-
+        .dvPlot(.datCmt, c("PRED", "IPRED"), TRUE) +
+        ggplot2::ggtitle(cmt, "log-scale DV vs PRED/IPRED")
+    } else if (.hasIpred) {
+      .lst[["dv_ipred_linear"]] <-
+        .dvPlot(.datCmt, "IPRED") +
+        ggplot2::ggtitle(cmt, "DV vs IPRED")
 
-      if (.hasCwres) {
-        .p1 <- .dvPlot(.dat0, c("CPRED", "IPRED")) +
-          ggplot2::ggtitle(.cmt, "DV vs CPRED/IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
+      .lst[["dv_ipred_log"]] <-
+        .dvPlot(.datCmt, "IPRED", TRUE) +
+        ggplot2::ggtitle(cmt, "log-scale DV vs IPRED")
+    }
 
-        .p1 <- .dvPlot(.dat0, c("CPRED", "IPRED"), TRUE) +
-          ggplot2::ggtitle(.cmt, "log-scale DV vs CPRED/IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
-      }
+    if (.hasCwres) {
+      .lst[["dv_cpred_linear"]] <-
+        .dvPlot(.datCmt, c("CPRED", "IPRED")) +
+        ggplot2::ggtitle(cmt, "DV vs CPRED/IPRED")
+
+      .lst[["dv_cpred_log"]] <-
+        .dvPlot(.datCmt, c("CPRED", "IPRED"), TRUE) +
+        ggplot2::ggtitle(cmt, "log-scale DV vs CPRED/IPRED")
+    }
 
-      if (.hasNpde) {
-        .p1 <- .dvPlot(.dat0, c("EPRED", "IPRED")) +
-          ggplot2::ggtitle(.cmt, "DV vs EPRED/IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
+    if (.hasNpde) {
+      .lst[["dv_epred_linear"]] <-
+        .dvPlot(.datCmt, c("EPRED", "IPRED")) +
+        ggplot2::ggtitle(cmt, "DV vs EPRED/IPRED")
 
-        .p1 <- .dvPlot(.dat0, c("EPRED", "IPRED"), TRUE) +
-          ggplot2::ggtitle(.cmt, "log-scale DV vs EPRED/IPRED")
-        .lst[[length(.lst) + 1]] <- .p1
-      }
+      .lst[["dv_epred_log"]] <-
+        .dvPlot(.datCmt, c("EPRED", "IPRED"), TRUE) +
+        ggplot2::ggtitle(cmt, "log-scale DV vs EPRED/IPRED")
+    }
 
-      for (x in c("IPRED", "PRED", "CPRED", "EPRED", "TIME", "tad")) {
-        if (any(names(.dat0) == x)) {
-          for (y in c("IWRES", "IRES", "RES", "CWRES", "NPD")) {
-            if (any(names(.dat0) == y)) {
-              if (y == "CWRES" && x %in% c("TIME", "CPRED")) {
-                .doIt <- TRUE
-              } else if (y == "NPD" && x %in% c("TIME", "EPRED")) {
-                .doIt <- TRUE
-              } else if (!(y %in% c("CWRES", "NPD"))) {
-                .doIt <- TRUE
-              }
-              if (.doIt) {
-                .p2 <- .scatterPlot(.dat0, c(x, y), .cmt, log = FALSE)
-                .lst[[length(.lst) + 1]] <- .p2
-                .p2 <- .scatterPlot(.dat0, c(x, y), .cmt, log = TRUE)
-                .lst[[length(.lst) + 1]] <- .p2
-              }
+    for (x in c("IPRED", "PRED", "CPRED", "EPRED", "TIME", "tad")) {
+      if (any(names(.datCmt) == x)) {
+        for (y in c("IWRES", "IRES", "RES", "CWRES", "NPD")) {
+          if (any(names(.datCmt) == y)) {
+            if (y == "CWRES" && x %in% c("TIME", "CPRED")) {
+              .doIt <- TRUE
+            } else if (y == "NPD" && x %in% c("TIME", "EPRED")) {
+              .doIt <- TRUE
+            } else if (!(y %in% c("CWRES", "NPD"))) {
+              .doIt <- TRUE
+            }
+            if (.doIt) {
+              .lst[[paste(y, x, "linear", sep = "_")]] <-
+                .scatterPlot(.datCmt, c(x, y), cmt, log = FALSE)
+              .lst[[paste(y, x, "log", sep = "_")]] <-
+                .scatterPlot(.datCmt, c(x, y), cmt, log = TRUE)
             }
           }
         }
       }
-      # .idPlot <- try(plot.nlmixr2AugPred(nlmixr2AugPred(object)));
-      # if (inherits(.idPlot, "try-error")){
-      .ids <- unique(.dat0$ID)
-      .s <- seq(1, length(.ids), by = 16)
-      .j <- 0
-      for (i in .s) {
-        .j <- .j + 1
-        .tmp <- .ids[seq(i, i + 15)]
-        .tmp <- .tmp[!is.na(.tmp)]
-        .d1 <- .dat0[.dat0$ID %in% .tmp, ]
+    }
+    # .idPlot <- try(plot.nlmixr2AugPred(nlmixr2AugPred(object)));
+    # if (inherits(.idPlot, "try-error")){
+    .ids <- unique(.datCmt$ID)
+    .s <- seq(1, length(.ids), by = 16)
+    .j <- 0
+    for (i in .s) {
+      .j <- .j + 1
+      .tmp <- .ids[seq(i, i + 15)]
+      .tmp <- .tmp[!is.na(.tmp)]
+      .d1 <- .datCmt[.datCmt$ID %in% .tmp, ]
 
-        .p3 <- ggplot2::ggplot(.d1, ggplot2::aes(x = .data$TIME, y = .data$DV)) +
-          ggplot2::geom_point() +
-          ggplot2::geom_line(ggplot2::aes(x = .data$TIME, y = .data$IPRED), col = "red", size = 1.2)
-        if (any(names(.d1) == "PRED")) {
-          .p3 <- .p3 + ggplot2::geom_line(ggplot2::aes(x = .data$TIME, y = .data$PRED), col = "blue", size = 1.2)
-        }
-        .p3 <- .p3 + ggplot2::facet_wrap(~ID) +
-          ggplot2::ggtitle(.cmt, sprintf("Individual Plots (%s of %s)", .j, length(.s))) +
-          rxode2::rxTheme()
-        if (any(names(.d1) == "lowerLim")) {
-          .p3 <-
-            .p3 +
-            geom_cens(ggplot2::aes(lower = .data$lowerLim, upper = .data$upperLim), fill = "purple")
-        }
-        .lst[[length(.lst) + 1]] <- .p3
+      .pIndividual <- ggplot2::ggplot(.d1, ggplot2::aes(x = .data$TIME, y = .data$DV)) +
+        ggplot2::geom_point() +
+        ggplot2::geom_line(ggplot2::aes(x = .data$TIME, y = .data$IPRED), col = "red", size = 1.2)
+      if (any(names(.d1) == "PRED")) {
+        .pIndividual <- .pIndividual + ggplot2::geom_line(ggplot2::aes(x = .data$TIME, y = .data$PRED), col = "blue", size = 1.2)
       }
-      .dat0$id2 <- factor(paste0("id: ", .dat0$ID, "; dose#: ", .dat0$dosenum))
-      .ids <- unique(.dat0$id2)
-      .s <- seq(1, length(.ids), by = 16)
-      .j <- 0
-      # for (i in .s) {
-      #   .j <- .j + 1
-      #   .tmp <- .ids[seq(i, i + 15)]
-      #   .tmp <- .tmp[!is.na(.tmp)]
-      #   .d1 <- .dat0[.dat0$id2 %in% .tmp, ]
-
-      #   .p3 <- ggplot2::ggplot(.d1, ggplot2::aes(x = tad, y = DV)) +
-      #     ggplot2::geom_point() +
-      #     ggplot2::geom_line(ggplot2::aes(x = tad, y = IPRED), col = "red", size = 1.2) +
-      #     ggplot2::geom_line(ggplot2::aes(x = tad, y = PRED), col = "blue", size = 1.2) +
-      #     ggplot2::facet_wrap(~id2) +
-      #     ggplot2::ggtitle(.cmt, sprintf("Individual TAD Plots (%s of %s)", .j, length(.s))) +
-      #     rxode2::rxTheme()
-      #   if (any(names(.d1) == "lowerLim")) {
-      #     .p3 <- .p3 + geom_cens(ggplot2::aes(lower=lowerLim, upper=upperLim), fill="purple")
-      #   }
-      #   .lst[[length(.lst) + 1]] <- .p3
-      # }
+      .pIndividual <- .pIndividual + ggplot2::facet_wrap(~ID) +
+        ggplot2::ggtitle(cmt, sprintf("Individual Plots (%s of %s)", .j, length(.s))) +
+        rxode2::rxTheme()
+      if (any(names(.d1) == "lowerLim")) {
+        .pIndividual <-
+          .pIndividual +
+          geom_cens(ggplot2::aes(lower = .data$lowerLim, upper = .data$upperLim), fill = "purple")
+      }
+      .lst[[paste("individual", i, sep = "_")]] <- .pIndividual
     }
+    .datCmt$id2 <- factor(paste0("id: ", .datCmt$ID, "; dose#: ", .datCmt$dosenum))
+    .ids <- unique(.datCmt$id2)
+    .s <- seq(1, length(.ids), by = 16)
+    .j <- 0
   }
-
-  # .id <- unique(.dat0$id2)
-  # if (grDevices::dev.cur() != 1){
-  #     .x  <- .lst
-  #     for (.i in seq_along(.x)){
-  #         plot(.x[[.i]])
-  #     }
-  # }
   class(.lst) <- "nlmixr2PlotList"
-  return(.lst)
+  .lst
 }
 
 #' @export
@@ -325,10 +309,9 @@ plot.nlmixr2FitCore <- function(x, ...) {
   stop("This is not a nlmixr2 data frame and cannot be plotted")
 }
 
-##' @export
+#' @export
 plot.nlmixr2FitCoreSilent <- plot.nlmixr2FitCore
 
-#'
 #' @title Produce trace-plot for fit if applicable
 #'
 #' @param x fit object
@@ -396,5 +379,5 @@ traceplot.nlmixr2FitCore <- function(x, ...) {
   }
 }
 
-##' @export
+#' @export
 traceplot.nlmixr2FitCoreSilent <- traceplot.nlmixr2FitCore
diff --git a/man/plot.nlmixr2FitData.Rd b/man/plot.nlmixr2FitData.Rd