From 8d21c9966305627f9ea8599503fb3709ca916f9b Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 13 Nov 2020 10:27:00 +0100
Subject: Fixed plotting of non-analytical solutions in saem

---
 R/plot.mixed.mmkin.R   |  8 +++-----
 R/saem.R               | 31 ++++++++++++++++++++-----------
 R/summary.saem.mmkin.R |  9 +++------
 3 files changed, 26 insertions(+), 22 deletions(-)

(limited to 'R')

diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 0d6b15c0..d8a1c2ac 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -75,8 +75,8 @@ plot.mixed.mmkin <- function(x,
     rownames(degparms_optim) <- ds_names
     degparms_optim_names <- setdiff(names(fit_1$par), names(fit_1$errparms))
     colnames(degparms_optim) <- degparms_optim_names
-    residual_type = ifelse(standardized, "iwres", "ires")
-    residuals <- x$so@results@predictions[[residual_type]]
+    residual_type = ifelse(standardized, "residual", "standardized")
+    residuals <- x$data[[residual_type]]
     degparms_pop <- x$so@results@fixed.effects
     names(degparms_pop) <- degparms_optim_names
   }
@@ -95,9 +95,7 @@ plot.mixed.mmkin <- function(x,
   odeini_names <- grep("_0$", degparms_all_names, value = TRUE)
   odeparms_names <- setdiff(degparms_all_names, odeini_names)
 
-  residual_type = ifelse(standardized, "iwres", "ires")
-
-  observed <- cbind(x$data,
+  observed <- cbind(x$data[c("ds", "name", "time", "value")],
     residual = residuals)
 
   solution_type = fit_1$solution_type
diff --git a/R/saem.R b/R/saem.R
index 377da903..18152baf 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -63,7 +63,7 @@
 #' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
 #'   A1 = mkinsub("SFO"))
 #' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO,
-#' # and compiled ODEs for FOMC, both are fast
+#' # and compiled ODEs for FOMC that are much slower
 #' f_mmkin <- mmkin(list(
 #'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
 #'   ds, quiet = TRUE)
@@ -75,11 +75,12 @@
 #' # We can use print, plot and summary methods to check the results
 #' print(f_saem_dfop_sfo)
 #' plot(f_saem_dfop_sfo)
-#' summary(f_saem_dfop_sfo, data = FALSE)
+#' summary(f_saem_dfop_sfo, data = TRUE)
 #'
 #' # Using a single core, the following takes about 6 minutes as we do not have an
 #' # analytical solution. Using 10 cores it is slower instead of faster
-#' #f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
+#' f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1)
+#' plot(f_saem_fomc)
 #' }
 #' @export
 saem <- function(object, control, ...) UseMethod("saem")
@@ -103,12 +104,6 @@ saem.mmkin <- function(object,
   }
   fit_time <- system.time({
     utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
-    f_pred <- try(saemix::saemix.predict(f_saemix), silent = TRUE)
-    if (!inherits(f_pred, "try-error")) {
-      f_saemix <- f_pred
-    } else {
-      warning("Creating predictions from the saemix model failed")
-    }
   })
   if (suppressPlot) {
     grDevices::dev.off()
@@ -121,6 +116,20 @@ saem.mmkin <- function(object,
     object[[1]]$transform_rates,
     object[[1]]$transform_fractions)
 
+  return_data <- nlme_data(object)
+
+  return_data$predicted <- f_saemix@model@model(
+    psi = saemix::psi(f_saemix),
+    id = as.numeric(return_data$ds),
+    xidep = return_data[c("time", "name")])
+
+  return_data <- transform(return_data,
+    residual = predicted - value,
+    std = sigma_twocomp(predicted,
+      f_saemix@results@respar[1], f_saemix@results@respar[2]))
+  return_data <- transform(return_data,
+    standardized = residual / std)
+
   result <- list(
     mkinmod = object[[1]]$mkinmod,
     mmkin = object,
@@ -132,7 +141,7 @@ saem.mmkin <- function(object,
     mean_dp_start = attr(m_saemix, "mean_dp_start"),
     bparms.optim = bparms_optim,
     bparms.fixed = object[[1]]$bparms.fixed,
-    data = nlme_data(object),
+    data = return_data,
     err_mod = object[[1]]$err_mod,
     date.fit = date(),
     saemixversion = as.character(utils::packageVersion("saemix")),
@@ -323,7 +332,7 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
       uid <- unique(id)
 
       res_list <- parallel::mclapply(uid, function(i) {
-          transparms_optim <- psi[i, ]
+          transparms_optim <- as.numeric(psi[i, ])
           names(transparms_optim) <- names(degparms_optim)
 
           odeini_optim <- transparms_optim[odeini_optim_parm_names]
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 536693a3..97f9f2da 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -151,14 +151,11 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
   err_mod = object$mmkin_orig[[1]]$err_mod
 
   object$diffs <- object$mkinmod$diffs
-  object$print_data <- data
+  object$print_data <- data # boolean: Should we print the data?
   so_pred <- object$so@results@predictions
 
-  object$data[["observed"]] <- object$data[["value"]]
-  object$data[["value"]] <- NULL
-  object$data[["predicted"]] <- so_pred$ipred
-  object$data[["residual"]] <- so_pred$ires
-  object$data[["standardized"]] <- so_pred$iwres
+  names(object$data)[4] <- "observed" # rename value to observed
+
   object$verbose <- verbose
 
   object$fixed <- object$mmkin_orig[[1]]$fixed
-- 
cgit v1.2.1