From dd80ab5d64bc6b56587b0542dc95e965a3a25590 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Sat, 24 Oct 2020 00:32:00 +0200
Subject: Plot method for nlme.mmkin objects

Update docs
---
 R/nlme.mmkin.R      |   2 +
 R/plot.mmkin.R      |  18 ++--
 R/plot.nlme.mmkin.R | 259 +++++++++++++++++++++++++++++++++++++---------------
 3 files changed, 196 insertions(+), 83 deletions(-)

(limited to 'R')

diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index c8a99d59..d4720e79 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -196,6 +196,8 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
 
   val <- do.call("nlme.formula", thisCall)
   val$mmkin_orig <- model
+  val$data <- thisCall[["data"]]
+  val$mkinmod <- model[[1]]$mkinmod
   class(val) <- c("nlme.mmkin", "nlme", "lme")
   return(val)
 }
diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R
index ac66e542..0523e4d3 100644
--- a/R/plot.mmkin.R
+++ b/R/plot.mmkin.R
@@ -54,12 +54,16 @@
 #'
 #' @export
 plot.mmkin <- function(x, main = "auto", legends = 1,
-                       resplot = c("time", "errmod"),
-                       standardized = FALSE,
-                       show_errmin = TRUE,
-                       errmin_var = "All data", errmin_digits = 3,
-                       cex = 0.7, rel.height.middle = 0.9,
-                       ymax = "auto", ...) {
+  resplot = c("time", "errmod"),
+  standardized = FALSE,
+  show_errmin = TRUE,
+  errmin_var = "All data", errmin_digits = 3,
+  cex = 0.7, rel.height.middle = 0.9,
+  ymax = "auto", ...)
+{
+
+  oldpar <- par(no.readonly = TRUE)
+
   n.m <- nrow(x)
   n.d <- ncol(x)
 
@@ -82,8 +86,6 @@ plot.mmkin <- function(x, main = "auto", legends = 1,
                   datasets = rownames(x))
   }
 
-  oldpar <- par(no.readonly = TRUE)
-
   # Set relative plot heights, so the first and the last plot are the norm
   # and the middle plots (if n.fits >2) are smaller by rel.height.middle
   rel.heights <- if (n.fits > 2) c(1, rep(rel.height.middle, n.fits - 2), 1)
diff --git a/R/plot.nlme.mmkin.R b/R/plot.nlme.mmkin.R
index 0f3ad715..084099ac 100644
--- a/R/plot.nlme.mmkin.R
+++ b/R/plot.nlme.mmkin.R
@@ -1,115 +1,224 @@
+if(getRversion() >= '2.15.1') utils::globalVariables("ds")
+
 #' Plot a fitted nonlinear mixed model obtained via an mmkin row object
 #'
 #' @param x An object of class \code{\link{nlme.mmkin}}
 #' @param i A numeric index to select datasets for which to plot the nlme fit,
 #'   in case plots get too large
 #' @param main The main title placed on the outer margin of the plot.
+#' @inheritParams plot.mkinfit
 #' @param legends An index for the fits for which legends should be shown.
-#' @param resplot Should the residuals plotted against time, using
-#'   \code{\link{mkinresplot}}, or as squared residuals against predicted
-#'   values, with the error model, using \code{\link{mkinerrplot}}.
-#' @param standardized Should the residuals be standardized? This option
-#'   is passed to \code{\link{mkinresplot}}, it only takes effect if
+#' @param standardized Should the residuals be standardized? Only takes effect if
 #'   `resplot = "time"`.
-#' @param show_errmin Should the chi2 error level be shown on top of the plots
-#'   to the left?
-#' @param errmin_var The variable for which the FOCUS chi2 error value should
-#'   be shown.
-#' @param errmin_digits The number of significant digits for rounding the FOCUS
-#'   chi2 error percentage.
 #' @param cex Passed to the plot functions and \code{\link{mtext}}.
 #' @param rel.height.middle The relative height of the middle plot, if more
 #'   than two rows of plots are shown.
-#' @param ymax Maximum y axis value for \code{\link{plot.mkinfit}}.
+#' @param ymax Vector of maximum y axis values
 #' @param \dots Further arguments passed to \code{\link{plot.mkinfit}} and
 #'   \code{\link{mkinresplot}}.
+#' @param resplot Should the residuals plotted against time or against
+#'   predicted values?
+#' @param col_ds Colors used for plotting the observed data and the
+#'   corresponding model prediction lines for the different datasets.
+#' @param pch_ds Symbols to be used for plotting the data.
+#' @param lty_ds Line types to be used for the model predictions.
 #' @importFrom stats coefficients
 #' @return The function is called for its side effect.
 #' @author Johannes Ranke
 #' @examples
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
-#'  function(x) subset(x$data[c("name", "time", "value")], name == "parent"))
-#' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1)
-#' #plot(f) # too many panels for pkgdown
-#' plot(f[, 3:4])
+#'  function(x) x$data[c("name", "time", "value")])
+#' names(ds) <- paste0("ds ", 6:10)
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+#'   A1 = mkinsub("SFO"), quiet = TRUE)
+#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, cores = 1)
+#' plot(f[, 3:4], standardized = TRUE)
 #' library(nlme)
-#' f_nlme <- nlme(f)
-#'
-#' #plot(f_nlme) # too many panels for pkgdown
-#' plot(f_nlme, 3:4)
+#' # For this fit we need to increase pnlsMaxiter, and we increase the
+#' # tolerance in order to speed up the fit for this example evaluation
+#' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-4))
+#' plot(f_nlme)
 #' @export
 plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig),
-  main = "auto", legends = 1,
-  resplot = c("time", "errmod"),
-  standardized = FALSE,
-  show_errmin = TRUE,
-  errmin_var = "All data", errmin_digits = 3,
+  main = rownames(x$mmkin_orig),
+  obs_vars = names(x$mkinmod$map),
+  standardized = TRUE,
+  xlab = "Time", ylab = "Observed",
+  xlim = range(x$data$time),
+  legends = 1,
+  lpos = "topright", inset = c(0.05, 0.05),
+  resplot = c("predicted", "time"),
+  ymax = "auto", maxabs = "auto",
   cex = 0.7, rel.height.middle = 0.9,
-  ymax = "auto", ...)
+  pch_ds = 1:length(i),
+  col_ds = pch_ds + 1,
+  lty_ds = col_ds,
+  frame = TRUE, ...)
 {
 
-  degparms_optim_nlme <- coefficients(x)
-  degparms_optim_names <- names(degparms_optim_nlme)
+  oldpar <- par(no.readonly = TRUE)
 
-  odeini_optim_names <- grep("_0$", degparms_optim_names, value = TRUE)
-  odeparms_optim_names <- setdiff(degparms_optim_names, odeini_optim_names)
+  fit_1 = x$mmkin_orig[[1]]
+  ds_names <- colnames(x$mmkin_orig)
 
-  fit_1 <- x$mmkin_orig[[1]]
+  degparms_optim <- coefficients(x)
+  degparms_optim_names <- names(degparms_optim)
+  degparms_fixed <- fit_1$fixed$value
+  names(degparms_fixed) <- rownames(fit_1$fixed)
+  degparms_all <- cbind(as.matrix(degparms_optim), fit_1$bparms.fixed)
+  degparms_all_names <- c(degparms_optim_names, names(degparms_fixed))
+  colnames(degparms_all) <- degparms_all_names
 
-  mkinfit_call <- as.list(fit_1$call)[-1]
-  mkinfit_call[["mkinmod"]] <- fit_1$mkinmod
+  odeini_names <- grep("_0$", degparms_all_names, value = TRUE)
+  odeparms_names <- setdiff(degparms_all_names, odeini_names)
 
-  ds <- lapply(x$mmkin_orig, function(x) {
-    data.frame(name = x$data$variable,
-      time = x$data$time,
-      value = x$data$observed)
-  })
+  residual_type = ifelse(standardized, "pearson", "response")
 
-  # This takes quite some time. This could be greatly reduced
-  # if the plot.mkinfit code would be imported and adapted,
-  # allowing also to overly plots of mmkin fits and nlme fits
-  mmkin_nlme <- lapply(i, function(a) {
+  observed <- cbind(x$data,
+    residual = residuals(x, type = residual_type))
 
-    degparms_optim <- as.numeric(degparms_optim_nlme[a, ])
-    names(degparms_optim) <- degparms_optim_names
+  n_plot_rows = length(obs_vars)
+  n_plots = n_plot_rows * 2
 
-    odeini_optim <- degparms_optim[odeini_optim_names]
-    names(odeini_optim) <- gsub("_0$", "", names(odeini_optim))
-
-    odeparms_optim_trans <- degparms_optim[odeparms_optim_names]
-    odeparms_optim <- backtransform_odeparms(odeparms_optim_trans,
-      fit_1$mkinmod,
-      transform_rates = fit_1$transform_rates,
-      transform_fractions = fit_1$transform_fractions)
+  # Set relative plot heights, so the first and the last plot are the norm
+  # and the middle plots (if n_plot_rows >2) are smaller by rel.height.middle
+  rel.heights <- if (n_plot_rows > 2) {
+    c(1, rep(rel.height.middle, n_plot_rows - 2), 1)
+  } else {
+    rep(1, n_plot_rows)
+  }
 
-    fit_a <- x$mmkin_orig[[a]]
+  layout_matrix = matrix(1:n_plots,
+    n_plot_rows, 2, byrow = TRUE)
+  layout(layout_matrix, heights = rel.heights)
 
-    state_ini <- fit_a$bparms.state
-    state_ini[names(odeini_optim)] <- odeini_optim
+  solution_type = fit_1$solution_type
 
-    odeparms <- fit_a$bparms.ode
-    odeparms[names(odeparms_optim)] <- odeparms_optim
+  outtimes <- sort(unique(c(x$data$time,
+    seq(xlim[1], xlim[2], length.out = 50))))
 
-    mkinfit_call[["observed"]] <- ds[[a]]
-    mkinfit_call[["parms.ini"]] <- odeparms
-    mkinfit_call[["state.ini"]] <- state_ini
+  pred_ds <- purrr::map_dfr(i, function(ds_i)   {
+    odeparms_trans <- degparms_all[ds_i, odeparms_names]
+    odeparms <- backtransform_odeparms(odeparms_trans,
+      x$mkinmod,
+      transform_rates = fit_1$transform_rates,
+      transform_fractions = fit_1$transform_fractions)
 
-    mkinfit_call[["control"]] <- list(iter.max = 0)
-    mkinfit_call[["quiet"]] <- TRUE
+    odeini <- degparms_all[ds_i, odeini_names]
+    names(odeini) <- gsub("_0", "", names(odeini))
 
-    res <- suppressWarnings(do.call("mkinfit", mkinfit_call))
-    return(res)
+    out <- mkinpredict(x$mkinmod, odeparms, odeini,
+      outtimes, solution_type = solution_type,
+      atol = fit_1$atol, rtol = fit_1$rtol)
+    return(cbind(as.data.frame(out), ds = ds_names[ds_i]))
   })
 
-  # Set dimensions with names and the class (mmkin)
-  attributes(mmkin_nlme) <- attributes(x$mmkin_orig[, i])
-
-  plot(mmkin_nlme, main = main, legends = legends,
-    resplot = resplot, standardized = standardized,
-    show_errmin = show_errmin,
-    errmin_var = errmin_var, errmin_digits = errmin_digits,
-    cex = cex,
-    rel.height.middle = rel.height.middle,
-    ymax = ymax, ...)
-
+  degparms_all_pop <- c(fixef(x), degparms_fixed)
+
+  odeparms_pop_trans <- degparms_all_pop[odeparms_names]
+  odeparms_pop <- backtransform_odeparms(odeparms_pop_trans,
+    x$mkinmod,
+    transform_rates = fit_1$transform_rates,
+    transform_fractions = fit_1$transform_fractions)
+
+  odeini_pop <- degparms_all_pop[odeini_names]
+  names(odeini_pop) <- gsub("_0", "", names(odeini_pop))
+
+  pred_pop <- as.data.frame(
+    mkinpredict(x$mkinmod, odeparms_pop, odeini_pop,
+      outtimes, solution_type = solution_type,
+      atol = fit_1$atol, rtol = fit_1$rtol))
+
+  resplot <- match.arg(resplot)
+
+  # Loop plot rows
+  for (plot_row in 1:n_plot_rows) {
+
+    obs_var <- obs_vars[plot_row]
+    observed_row <- subset(observed, name == obs_var)
+
+    # Set ylim to sensible default, or use ymax
+    if (identical(ymax, "auto")) {
+      ylim_row = c(0,
+        max(c(observed_row$value, pred_ds[[obs_var]]), na.rm = TRUE))
+    } else {
+      ylim_row = c(0, ymax[plot_row])
+    }
+
+    # Margins for top row of plots when we have more than one row
+    # Reduce bottom margin by 2.1 - hides x axis legend
+    if (plot_row == 1 & n_plot_rows > 1) {
+      par(mar = c(3.0, 4.1, 4.1, 2.1))
+    }
+
+    # Margins for middle rows of plots, if any
+    if (plot_row > 1 & plot_row < n_plot_rows) {
+      # Reduce top margin by 2 after the first plot as we have no main title,
+      # reduced plot height, therefore we need rel.height.middle in the layout
+      par(mar = c(3.0, 4.1, 2.1, 2.1))
+    }
+
+    # Margins for bottom row of plots when we have more than one row
+    if (plot_row == n_plot_rows & n_plot_rows > 1) {
+      # Restore bottom margin for last plot to show x axis legend
+      par(mar = c(5.1, 4.1, 2.1, 2.1))
+    }
+
+    plot(pred_pop$time, pred_pop[[obs_var]],
+      main = obs_var,
+      type = "l", lwd = 2,
+      xlim = xlim, ylim = ylim_row,
+      xlab = xlab, ylab = ylab, frame = frame,
+      cex = cex)
+
+    for (ds_i in seq_along(i)) {
+      points(subset(observed_row, ds == ds_names[ds_i], c("time", "value")),
+        col = col_ds[ds_i], pch = pch_ds[ds_i])
+      lines(subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)),
+        col = col_ds[ds_i], lty = lty_ds[ds_i])
+    }
+
+    if (plot_row %in% legends) {
+      legend(lpos, inset = inset,
+        legend = c("Population mean", ds_names[i]),
+        lty = c(1, lty_ds), lwd = c(2, rep(1, length(i))),
+        col = c(1, col_ds),
+        pch = c(NA, pch_ds))
+    }
+
+    if (identical(maxabs, "auto")) {
+      maxabs = max(abs(observed_row$residual), na.rm = TRUE)
+    }
+
+
+    if (identical(resplot, "time")) {
+      plot(0, type = "n", xlim = xlim, xlab = "Time",
+        main = obs_var,
+        ylim = c(-1.2 * maxabs, 1.2 * maxabs),
+        ylab = if (standardized) "Standardized residual" else "Residual")
+      abline(h = 0, lty = 2)
+
+      for (ds_i in seq_along(i)) {
+        points(subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
+          col = col_ds[ds_i], pch = pch_ds[ds_i])
+      }
+    }
+
+    if (identical(resplot, "predicted")) {
+      plot(0, type = "n",
+        main = obs_var,
+        xlim = c(0, max(pred_ds[[obs_var]])),
+        xlab = "Predicted",
+        ylim = c(-1.2 * maxabs, 1.2 * maxabs),
+        ylab = if (standardized) "Standardized residual" else "Residual")
+
+      for (ds_i in seq_along(i)) {
+        observed_row_ds <- merge(
+          subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
+          subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)))
+        points(observed_row_ds[c(3, 2)],
+          col = col_ds[ds_i], pch = pch_ds[ds_i])
+      }
+    }
+  }
 }
-- 
cgit v1.2.1