From 48c463680b51fa767b4cd7bd62865f192d0354ac Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Sat, 6 Feb 2021 18:30:32 +0100
Subject: Reintroduce interface to saemix

Also after the upgrade from buster to bullseye of my local system, some
test results for saemix have changed.
---
 R/endpoints.R          |   8 +-
 R/plot.mixed.mmkin.R   |  23 ++-
 R/saem.R               | 512 +++++++++++++++++++++++++++++++++++++++++++++++++
 R/summary.saem.mmkin.R | 268 ++++++++++++++++++++++++++
 4 files changed, 806 insertions(+), 5 deletions(-)
 create mode 100644 R/saem.R
 create mode 100644 R/summary.saem.mmkin.R

(limited to 'R')

diff --git a/R/endpoints.R b/R/endpoints.R
index b5872e68..f1f47581 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -10,8 +10,8 @@
 #' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
 #' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
 #'
-#' @param fit An object of class [mkinfit] or [nlme.mmkin] 
-#'  or another object that has list components
+#' @param fit An object of class [mkinfit], [nlme.mmkin] or
+#'  [saem.mmkin]. Or another object that has list components
 #'  mkinmod containing an [mkinmod] degradation model, and two numeric vectors,
 #'  bparms.optim and bparms.fixed, that contain parameter values
 #'  for that model.
@@ -20,8 +20,8 @@
 #'   and, if applicable, a vector of formation fractions named ff
 #'   and, if the SFORB model was in use, a vector of eigenvalues
 #'   of these SFORB models, equivalent to DFOP rate constants
-#' @note The function is used internally by [summary.mkinfit]
-#'   and [summary.nlme.mmkin]
+#' @note The function is used internally by [summary.mkinfit],
+#'   [summary.nlme.mmkin] and [summary.saem.mmkin].
 #' @author Johannes Ranke
 #' @examples
 #'
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 5a0b7412..1674d855 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -2,7 +2,7 @@ utils::globalVariables("ds")
 
 #' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
 #'
-#' @param x An object of class [mixed.mmkin], [nlme.mmkin]
+#' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin]
 #' @param i A numeric index to select datasets for which to plot the individual predictions,
 #'   in case plots get too large
 #' @inheritParams plot.mkinfit
@@ -39,6 +39,15 @@ utils::globalVariables("ds")
 #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
 #' plot(f_nlme)
 #'
+#' f_saem <- saem(f, transformations = "saemix")
+#' plot(f_saem)
+#'
+#' # We can overlay the two variants if we generate predictions
+#' pred_nlme <- mkinpredict(dfop_sfo,
+#'   f_nlme$bparms.optim[-1],
+#'   c(parent = f_nlme$bparms.optim[[1]], A1 = 0),
+#'   seq(0, 180, by = 0.2))
+#' plot(f_saem, pred_over = list(nlme = pred_nlme))
 #' }
 #' @export
 plot.mixed.mmkin <- function(x,
@@ -82,6 +91,18 @@ plot.mixed.mmkin <- function(x,
       type = ifelse(standardized, "pearson", "response"))
   }
 
+  if (inherits(x, "saem.mmkin")) {
+    if (x$transformations == "saemix") backtransform = FALSE
+    degparms_i <- saemix::psi(x$so)
+    rownames(degparms_i) <- ds_names
+    degparms_i_names <- setdiff(x$so@results@name.fixed, names(fit_1$errparms))
+    colnames(degparms_i) <- degparms_i_names
+    residual_type = ifelse(standardized, "standardized", "residual")
+    residuals <- x$data[[residual_type]]
+    degparms_pop <- x$so@results@fixed.effects
+    names(degparms_pop) <- degparms_i_names
+  }
+
   degparms_fixed <- fit_1$fixed$value
   names(degparms_fixed) <- rownames(fit_1$fixed)
   degparms_all <- cbind(as.matrix(degparms_i),
diff --git a/R/saem.R b/R/saem.R
new file mode 100644
index 00000000..fd2a77b4
--- /dev/null
+++ b/R/saem.R
@@ -0,0 +1,512 @@
+utils::globalVariables(c("predicted", "std"))
+
+#' Fit nonlinear mixed models with SAEM
+#'
+#' This function uses [saemix::saemix()] as a backend for fitting nonlinear mixed
+#' effects models created from [mmkin] row objects using the Stochastic Approximation
+#' Expectation Maximisation algorithm (SAEM).
+#'
+#' An mmkin row object is essentially a list of mkinfit objects that have been
+#' obtained by fitting the same model to a list of datasets using [mkinfit].
+#'
+#' Starting values for the fixed effects (population mean parameters, argument
+#' psi0 of [saemix::saemixModel()] are the mean values of the parameters found
+#' using [mmkin].
+#'
+#' @param object An [mmkin] row object containing several fits of the same
+#'   [mkinmod] model to different datasets
+#' @param verbose Should we print information about created objects of
+#'   type [saemix::SaemixModel] and [saemix::SaemixData]?
+#' @param transformations Per default, all parameter transformations are done
+#'   in mkin. If this argument is set to 'saemix', parameter transformations
+#'   are done in 'saemix' for the supported cases. Currently this is only
+#'   supported in cases where the initial concentration of the parent is not fixed,
+#'   SFO or DFOP is used for the parent and there is either no metabolite or one.
+#' @param degparms_start Parameter values given as a named numeric vector will
+#'   be used to override the starting values obtained from the 'mmkin' object.
+#' @param solution_type Possibility to specify the solution type in case the
+#'   automatic choice is not desired
+#' @param quiet Should we suppress the messages saemix prints at the beginning
+#'   and the end of the optimisation process?
+#' @param control Passed to [saemix::saemix]
+#' @param \dots Further parameters passed to [saemix::saemixModel].
+#' @return An S3 object of class 'saem.mmkin', containing the fitted
+#'   [saemix::SaemixObject] as a list component named 'so'. The
+#'   object also inherits from 'mixed.mmkin'.
+#' @seealso [summary.saem.mmkin] [plot.mixed.mmkin]
+#' @examples
+#' \dontrun{
+#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
+#'  function(x) subset(x$data[c("name", "time", "value")]))
+#' names(ds) <- paste("Dataset", 6:10)
+#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds,
+#'   state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
+#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
+#'
+#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
+#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
+#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
+#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+#'
+#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
+#' # functions from saemix
+#' library(saemix)
+#' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so))
+#' plot(f_saem_fomc$so, plot.type = "convergence")
+#' plot(f_saem_fomc$so, plot.type = "individual.fit")
+#' plot(f_saem_fomc$so, plot.type = "npde")
+#' plot(f_saem_fomc$so, plot.type = "vpc")
+#'
+#' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
+#' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
+#' compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
+#'
+#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+#'   A1 = mkinsub("SFO"))
+#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+#'   A1 = mkinsub("SFO"))
+#' 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 that are much slower
+#' f_mmkin <- mmkin(list(
+#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#'   ds, quiet = TRUE)
+#' # saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds
+#' # each on this system, as we use analytical solutions written for saemix.
+#' # When using the analytical solutions written for mkin this took around
+#' # four minutes
+#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
+#' f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
+#' # 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 = TRUE)
+#'
+#' # The following takes about 6 minutes
+#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
+#' #  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+#'
+#' #saemix::compare.saemix(list(
+#' #  f_saem_dfop_sfo$so,
+#' #  f_saem_dfop_sfo_deSolve$so))
+#'
+#' # If the model supports it, we can also use eigenvalue based solutions, which
+#' # take a similar amount of time
+#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen",
+#' #  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+#' }
+#' @export
+saem <- function(object, ...) UseMethod("saem")
+
+#' @rdname saem
+#' @export
+saem.mmkin <- function(object,
+  transformations = c("mkin", "saemix"),
+  degparms_start = numeric(),
+  solution_type = "auto",
+  control = list(displayProgress = FALSE, print = FALSE,
+    save = FALSE, save.graphs = FALSE),
+  verbose = FALSE, quiet = FALSE, ...)
+{
+  transformations <- match.arg(transformations)
+  m_saemix <- saemix_model(object, verbose = verbose,
+    degparms_start = degparms_start, solution_type = solution_type,
+    transformations = transformations, ...)
+  d_saemix <- saemix_data(object, verbose = verbose)
+
+  fit_time <- system.time({
+    utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
+  })
+
+  transparms_optim <- f_saemix@results@fixed.effects
+  names(transparms_optim) <- f_saemix@results@name.fixed
+
+  if (transformations == "mkin") {
+    bparms_optim <- backtransform_odeparms(transparms_optim,
+      object[[1]]$mkinmod,
+      object[[1]]$transform_rates,
+      object[[1]]$transform_fractions)
+  } else {
+    bparms_optim <- transparms_optim
+  }
+
+  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,
+    solution_type = object[[1]]$solution_type,
+    transformations = transformations,
+    transform_rates = object[[1]]$transform_rates,
+    transform_fractions = object[[1]]$transform_fractions,
+    so = f_saemix,
+    time = fit_time,
+    mean_dp_start = attr(m_saemix, "mean_dp_start"),
+    bparms.optim = bparms_optim,
+    bparms.fixed = object[[1]]$bparms.fixed,
+    data = return_data,
+    err_mod = object[[1]]$err_mod,
+    date.fit = date(),
+    saemixversion = as.character(utils::packageVersion("saemix")),
+    mkinversion = as.character(utils::packageVersion("mkin")),
+    Rversion = paste(R.version$major, R.version$minor, sep=".")
+  )
+
+  class(result) <- c("saem.mmkin", "mixed.mmkin")
+  return(result)
+}
+
+#' @export
+#' @rdname saem
+#' @param x An saem.mmkin object to print
+#' @param digits Number of digits to use for printing
+print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+  cat( "Kinetic nonlinear mixed-effects model fit by SAEM" )
+  cat("\nStructural model:\n")
+  diffs <- x$mmkin[[1]]$mkinmod$diffs
+  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs)
+  writeLines(strwrap(nice_diffs, exdent = 11))
+  cat("\nData:\n")
+  cat(nrow(x$data), "observations of",
+    length(unique(x$data$name)), "variable(s) grouped in",
+    length(unique(x$data$ds)), "datasets\n")
+
+  cat("\nLikelihood computed by importance sampling\n")
+  print(data.frame(
+      AIC = AIC(x$so, type = "is"),
+      BIC = BIC(x$so, type = "is"),
+      logLik = logLik(x$so, type = "is"),
+      row.names = " "), digits = digits)
+
+  cat("\nFitted parameters:\n")
+  conf.int <- x$so@results@conf.int[c("estimate", "lower", "upper")]
+  rownames(conf.int) <- x$so@results@conf.int[["name"]]
+  print(conf.int, digits = digits)
+
+  invisible(x)
+}
+
+#' @rdname saem
+#' @return An [saemix::SaemixModel] object.
+#' @export
+saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
+  degparms_start = numeric(), verbose = FALSE, ...)
+{
+  if (nrow(object) > 1) stop("Only row objects allowed")
+
+  mkin_model <- object[[1]]$mkinmod
+
+  degparms_optim <-  mean_degparms(object)
+  if (transformations == "saemix") {
+    degparms_optim <- backtransform_odeparms(degparms_optim,
+      object[[1]]$mkinmod,
+      object[[1]]$transform_rates,
+      object[[1]]$transform_fractions)
+  }
+  degparms_fixed <- object[[1]]$bparms.fixed
+
+  # Transformations are done in the degradation function
+  transform.par = rep(0, length(degparms_optim))
+
+  odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
+  odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE)
+
+  odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names)
+  odeparms_fixed <- degparms_fixed[odeparms_fixed_names]
+
+  odeini_fixed <- degparms_fixed[odeini_fixed_parm_names]
+  names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names)
+
+  model_function <- FALSE
+
+  # Model functions with analytical solutions
+  # Fixed parameters, use_of_ff = "min" and turning off sinks currently not supported here
+  # In general, we need to consider exactly how the parameters in mkinfit were specified,
+  # as the parameters are currently mapped by position in these solutions
+  sinks <- sapply(mkin_model$spec, function(x) x$sink)
+  if (length(odeparms_fixed) == 0 & mkin_model$use_of_ff == "max" & all(sinks)) {
+    # Parent only
+    if (length(mkin_model$spec) == 1) {
+      parent_type <- mkin_model$spec[[1]]$type
+      if (length(odeini_fixed) == 1) {
+        if (parent_type == "SFO") {
+          stop("saemix needs at least two parameters to work on.")
+        }
+        if (parent_type == "FOMC") {
+          model_function <- function(psi, id, xidep) {
+            odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1])
+          }
+        }
+        if (parent_type == "DFOP") {
+          model_function <- function(psi, id, xidep) {
+            g <- plogis(psi[id, 3])
+            t <- xidep[, "time"]
+            odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) +
+              (1 - g) * exp(- exp(psi[id, 2]) * t))
+          }
+        }
+        if (parent_type == "HS") {
+          model_function <- function(psi, id, xidep) {
+            tb <- exp(psi[id, 3])
+            t <- xidep[, "time"]
+            k1 = exp(psi[id, 1])
+            odeini_fixed * ifelse(t <= tb,
+              exp(- k1 * t),
+              exp(- k1 * tb) * exp(- exp(psi[id, 2]) * (t - tb)))
+          }
+        }
+      } else {
+        if (parent_type == "SFO") {
+          if (transformations == "mkin") {
+            model_function <- function(psi, id, xidep) {
+              psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"])
+            }
+          } else {
+            model_function <- function(psi, id, xidep) {
+              psi[id, 1] * exp( - psi[id, 2] * xidep[, "time"])
+            }
+            transform.par = c(0, 1)
+          }
+        }
+        if (parent_type == "FOMC") {
+          model_function <- function(psi, id, xidep) {
+            psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2])
+          }
+        }
+        if (parent_type == "DFOP") {
+          if (transformations == "mkin") {
+            model_function <- function(psi, id, xidep) {
+              g <- plogis(psi[id, 4])
+              t <- xidep[, "time"]
+              psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) +
+                (1 - g) * exp(- exp(psi[id, 3]) * t))
+            }
+          } else {
+            model_function <- function(psi, id, xidep) {
+              g <- psi[id, 4]
+              t <- xidep[, "time"]
+              psi[id, 1] * (g * exp(- psi[id, 2] * t) +
+                (1 - g) * exp(- psi[id, 3] * t))
+            }
+            transform.par = c(0, 1, 1, 3)
+          }
+        }
+        if (parent_type == "HS") {
+          model_function <- function(psi, id, xidep) {
+            tb <- exp(psi[id, 4])
+            t <- xidep[, "time"]
+            k1 = exp(psi[id, 2])
+            psi[id, 1] * ifelse(t <= tb,
+              exp(- k1 * t),
+              exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb)))
+          }
+        }
+      }
+    }
+
+    # Parent with one metabolite
+    # Parameter names used in the model functions are as in
+    # https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb
+    types <- unname(sapply(mkin_model$spec, function(x) x$type))
+    if (length(mkin_model$spec) == 2 &! "SFORB" %in% types ) {
+      # Initial value for the metabolite (n20) must be fixed
+      if (names(odeini_fixed) == names(mkin_model$spec)[2]) {
+        n20 <- odeini_fixed
+        parent_name <- names(mkin_model$spec)[1]
+        if (identical(types, c("SFO", "SFO"))) {
+          if (transformations == "mkin") {
+            model_function <- function(psi, id, xidep) {
+              t <- xidep[, "time"]
+              n10 <- psi[id, 1]
+              k1 <- exp(psi[id, 2])
+              k2 <- exp(psi[id, 3])
+              f12 <- plogis(psi[id, 4])
+              ifelse(xidep[, "name"] == parent_name,
+                n10 * exp(- k1 * t),
+                (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
+                  (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
+              )
+            }
+          } else {
+            model_function <- function(psi, id, xidep) {
+              t <- xidep[, "time"]
+              n10 <- psi[id, 1]
+              k1 <- psi[id, 2]
+              k2 <- psi[id, 3]
+              f12 <- psi[id, 4]
+              ifelse(xidep[, "name"] == parent_name,
+                n10 * exp(- k1 * t),
+                (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
+                  (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
+              )
+            }
+            transform.par = c(0, 1, 1, 3)
+          }
+        }
+        if (identical(types, c("DFOP", "SFO"))) {
+          if (transformations == "mkin") {
+            model_function <- function(psi, id, xidep) {
+              t <- xidep[, "time"]
+              n10 <- psi[id, 1]
+              k2 <- exp(psi[id, 2])
+              f12 <- plogis(psi[id, 3])
+              l1 <- exp(psi[id, 4])
+              l2 <- exp(psi[id, 5])
+              g <- plogis(psi[id, 6])
+              ifelse(xidep[, "name"] == parent_name,
+                n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
+                ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
+                  (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
+                  ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
+                      ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
+                        f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
+                  ((l1 - k2) * l2 - k2 * l1 + k2^2)
+              )
+            }
+          } else {
+            model_function <- function(psi, id, xidep) {
+              t <- xidep[, "time"]
+              n10 <- psi[id, 1]
+              k2 <- psi[id, 2]
+              f12 <- psi[id, 3]
+              l1 <- psi[id, 4]
+              l2 <- psi[id, 5]
+              g <- psi[id, 6]
+              ifelse(xidep[, "name"] == parent_name,
+                n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
+                ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
+                  (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
+                  ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
+                      ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
+                        f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
+                  ((l1 - k2) * l2 - k2 * l1 + k2^2)
+              )
+            }
+            transform.par = c(0, 1, 3, 1, 1, 3)
+          }
+        }
+      }
+    }
+  }
+
+  if (is.function(model_function) & solution_type == "auto") {
+    solution_type = "analytical saemix"
+  } else {
+
+    if (solution_type == "auto")
+      solution_type <- object[[1]]$solution_type
+
+    model_function <- function(psi, id, xidep) {
+
+      uid <- unique(id)
+
+      res_list <- lapply(uid, function(i) {
+
+        transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict
+        names(transparms_optim) <- names(degparms_optim)
+
+        odeini_optim <- transparms_optim[odeini_optim_parm_names]
+        names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names)
+
+        odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)]
+
+        ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names)
+        odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model,
+          transform_rates = object[[1]]$transform_rates,
+          transform_fractions = object[[1]]$transform_fractions)
+        odeparms <- c(odeparms_optim, odeparms_fixed)
+
+        xidep_i <- subset(xidep, id == i)
+
+        if (solution_type == "analytical") {
+          out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms)
+        } else {
+
+          i_time <- xidep_i$time
+          i_name <- xidep_i$name
+
+          out_wide <- mkinpredict(mkin_model,
+            odeparms = odeparms, odeini = odeini,
+            solution_type = solution_type,
+            outtimes = sort(unique(i_time)),
+            na_stop = FALSE
+          )
+
+          out_index <- cbind(as.character(i_time), as.character(i_name))
+          out_values <- out_wide[out_index]
+        }
+        return(out_values)
+      })
+      res <- unlist(res_list)
+      return(res)
+    }
+  }
+
+  error.model <- switch(object[[1]]$err_mod,
+    const = "constant",
+    tc = "combined",
+    obs = "constant")
+
+  if (object[[1]]$err_mod == "obs") {
+    warning("The error model 'obs' (variance by variable) can currently not be transferred to an saemix model")
+  }
+
+  error.init <- switch(object[[1]]$err_mod,
+    const = c(a = mean(sapply(object, function(x) x$errparms)), b = 1),
+    tc = c(a = mean(sapply(object, function(x) x$errparms[1])),
+      b = mean(sapply(object, function(x) x$errparms[2]))),
+    obs = c(a = mean(sapply(object, function(x) x$errparms)), b = 1))
+
+  degparms_psi0 <- degparms_optim
+  degparms_psi0[names(degparms_start)] <- degparms_start
+  psi0_matrix <- matrix(degparms_psi0, nrow = 1)
+  colnames(psi0_matrix) <- names(degparms_psi0)
+
+  res <- saemix::saemixModel(model_function,
+    psi0 = psi0_matrix,
+    "Mixed model generated from mmkin object",
+    transform.par = transform.par,
+    error.model = error.model,
+    verbose = verbose
+  )
+  attr(res, "mean_dp_start") <- degparms_optim
+  return(res)
+}
+
+#' @rdname saem
+#' @return An [saemix::SaemixData] object.
+#' @export
+saemix_data <- function(object, verbose = FALSE, ...) {
+  if (nrow(object) > 1) stop("Only row objects allowed")
+  ds_names <- colnames(object)
+
+  ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
+  names(ds_list) <- ds_names
+  ds_saemix_all <- purrr::map_dfr(ds_list, function(x) x, .id = "ds")
+  ds_saemix <- data.frame(ds = ds_saemix_all$ds,
+    name = as.character(ds_saemix_all$variable),
+    time = ds_saemix_all$time,
+    value = ds_saemix_all$observed,
+    stringsAsFactors = FALSE)
+
+  res <- saemix::saemixData(ds_saemix,
+    name.group = "ds",
+    name.predictors = c("time", "name"),
+    name.response = "value",
+    verbose = verbose,
+    ...)
+  return(res)
+}
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
new file mode 100644
index 00000000..e92c561c
--- /dev/null
+++ b/R/summary.saem.mmkin.R
@@ -0,0 +1,268 @@
+#' Summary method for class "saem.mmkin"
+#'
+#' Lists model equations, initial parameter values, optimised parameters
+#' for fixed effects (population), random effects (deviations from the
+#' population mean) and residual error model, as well as the resulting
+#' endpoints such as formation fractions and DT50 values. Optionally
+#' (default is FALSE), the data are listed in full.
+#'
+#' @param object an object of class [saem.mmkin]
+#' @param x an object of class [summary.saem.mmkin]
+#' @param data logical, indicating whether the full data should be included in
+#'   the summary.
+#' @param verbose Should the summary be verbose?
+#' @param distimes logical, indicating whether DT50 and DT90 values should be
+#'   included.
+#' @param digits Number of digits to use for printing
+#' @param \dots optional arguments passed to methods like \code{print}.
+#' @return The summary function returns a list based on the [saemix::SaemixObject]
+#'   obtained in the fit, with at least the following additional components
+#'   \item{saemixversion, mkinversion, Rversion}{The saemix, mkin and R versions used}
+#'   \item{date.fit, date.summary}{The dates where the fit and the summary were
+#'     produced}
+#'   \item{diffs}{The differential equations used in the degradation model}
+#'   \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
+#'   \item{data}{The data}
+#'   \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals}
+#'   \item{confint_back}{Backtransformed parameters, with confidence intervals if available}
+#'   \item{confint_errmod}{Error model parameters with confidence intervals}
+#'   \item{ff}{The estimated formation fractions derived from the fitted
+#'      model.}
+#'   \item{distimes}{The DT50 and DT90 values for each observed variable.}
+#'   \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
+#'   The print method is called for its side effect, i.e. printing the summary.
+#' @importFrom stats predict vcov
+#' @author Johannes Ranke for the mkin specific parts
+#'   saemix authors for the parts inherited from saemix.
+#' @examples
+#' # Generate five datasets following DFOP-SFO kinetics
+#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"),
+#'  m1 = mkinsub("SFO"), quiet = TRUE)
+#' set.seed(1234)
+#' k1_in <- rlnorm(5, log(0.1), 0.3)
+#' k2_in <- rlnorm(5, log(0.02), 0.3)
+#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3))
+#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3))
+#' k_m1_in <- rlnorm(5, log(0.02), 0.3)
+#'
+#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) {
+#'   mkinpredict(dfop_sfo,
+#'     c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1),
+#'     c(parent = 100, m1 = 0),
+#'     sampling_times)
+#' }
+#'
+#' ds_mean_dfop_sfo <- lapply(1:5, function(i) {
+#'   mkinpredict(dfop_sfo,
+#'     c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i],
+#'       f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]),
+#'     c(parent = 100, m1 = 0),
+#'     sampling_times)
+#' })
+#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5)
+#'
+#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) {
+#'   add_err(ds,
+#'     sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
+#'     n = 1)[[1]]
+#' })
+#'
+#' \dontrun{
+#' # Evaluate using mmkin and saem
+#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
+#'   quiet = TRUE, error_model = "tc", cores = 5)
+#' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo)
+#' summary(f_saem_dfop_sfo, data = TRUE)
+#' }
+#'
+#' @export
+summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) {
+
+  mod_vars <- names(object$mkinmod$diffs)
+
+  pnames <- names(object$mean_dp_start)
+  np <- length(pnames)
+
+  conf.int <- object$so@results@conf.int
+  rownames(conf.int) <- conf.int$name
+  confint_trans <- as.matrix(conf.int[pnames, c("estimate", "lower", "upper")])
+  colnames(confint_trans)[1] <- "est."
+
+  # In case objects were produced by earlier versions of saem
+  if (is.null(object$transformations)) object$transformations <- "mkin"
+
+  if (object$transformations == "mkin") {
+    bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
+      object$transform_rates, object$transform_fractions)
+    bpnames <- names(bp)
+
+    # Transform boundaries of CI for one parameter at a time,
+    # with the exception of sets of formation fractions (single fractions are OK).
+    f_names_skip <- character(0)
+    for (box in mod_vars) { # Figure out sets of fractions to skip
+      f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
+      n_paths <- length(f_names)
+      if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
+    }
+
+    confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
+      dimnames = list(bpnames, colnames(confint_trans)))
+    confint_back[, "est."] <- bp
+
+    for (pname in pnames) {
+      if (!pname %in% f_names_skip) {
+        par.lower <- confint_trans[pname, "lower"]
+        par.upper <- confint_trans[pname, "upper"]
+        names(par.lower) <- names(par.upper) <- pname
+        bpl <- backtransform_odeparms(par.lower, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        bpu <- backtransform_odeparms(par.upper, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        confint_back[names(bpl), "lower"] <- bpl
+        confint_back[names(bpu), "upper"] <- bpu
+      }
+    }
+  } else {
+    confint_back <- confint_trans
+  }
+
+  #  Correlation of fixed effects (inspired by summary.nlme)
+  varFix <- vcov(object$so)[1:np, 1:np]
+  stdFix <- sqrt(diag(varFix))
+  object$corFixed <- array(
+    t(varFix/stdFix)/stdFix,
+    dim(varFix),
+    list(pnames, pnames))
+
+  # Random effects
+  rnames <- paste0("SD.", pnames)
+  confint_ranef <- as.matrix(conf.int[rnames, c("estimate", "lower", "upper")])
+  colnames(confint_ranef)[1] <- "est."
+
+  # Error model
+  enames <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1")
+  confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")])
+  colnames(confint_errmod)[1] <- "est."
+
+
+  object$confint_trans <- confint_trans
+  object$confint_ranef <- confint_ranef
+  object$confint_errmod <- confint_errmod
+  object$confint_back <- confint_back
+
+  object$date.summary = date()
+  object$use_of_ff = object$mkinmod$use_of_ff
+  object$error_model_algorithm = object$mmkin_orig[[1]]$error_model_algorithm
+  err_mod = object$mmkin_orig[[1]]$err_mod
+
+  object$diffs <- object$mkinmod$diffs
+  object$print_data <- data # boolean: Should we print the data?
+  so_pred <- object$so@results@predictions
+
+  names(object$data)[4] <- "observed" # rename value to observed
+
+  object$verbose <- verbose
+
+  object$fixed <- object$mmkin_orig[[1]]$fixed
+  object$AIC = AIC(object$so)
+  object$BIC = BIC(object$so)
+  object$logLik = logLik(object$so, method = "is")
+
+  ep <- endpoints(object)
+  if (length(ep$ff) != 0)
+    object$ff <- ep$ff
+  if (distimes) object$distimes <- ep$distimes
+  if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB
+  class(object) <- c("summary.saem.mmkin")
+  return(object)
+}
+
+#' @rdname summary.saem.mmkin
+#' @export
+print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
+  cat("saemix version used for fitting:     ", x$saemixversion, "\n")
+  cat("mkin version used for pre-fitting: ", x$mkinversion, "\n")
+  cat("R version used for fitting:        ", x$Rversion, "\n")
+
+  cat("Date of fit:    ", x$date.fit, "\n")
+  cat("Date of summary:", x$date.summary, "\n")
+
+  cat("\nEquations:\n")
+  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
+  writeLines(strwrap(nice_diffs, exdent = 11))
+
+  cat("\nData:\n")
+  cat(nrow(x$data), "observations of",
+    length(unique(x$data$name)), "variable(s) grouped in",
+    length(unique(x$data$ds)), "datasets\n")
+
+  cat("\nModel predictions using solution type", x$solution_type, "\n")
+
+  cat("\nFitted in", x$time[["elapsed"]],  "s using", paste(x$so@options$nbiter.saemix, collapse = ", "), "iterations\n")
+
+  cat("\nVariance model: ")
+  cat(switch(x$err_mod,
+    const = "Constant variance",
+    obs = "Variance unique to each observed variable",
+    tc = "Two-component variance function"), "\n")
+
+  cat("\nMean of starting values for individual parameters:\n")
+  print(x$mean_dp_start, digits = digits)
+
+  cat("\nFixed degradation parameter values:\n")
+  if(length(x$fixed$value) == 0) cat("None\n")
+  else print(x$fixed, digits = digits)
+
+  cat("\nResults:\n\n")
+  cat("Likelihood computed by importance sampling\n")
+  print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
+      row.names = " "), digits = digits)
+
+  cat("\nOptimised parameters:\n")
+  print(x$confint_trans, digits = digits)
+
+  if (nrow(x$confint_trans) > 1) {
+    corr <- x$corFixed
+    class(corr) <- "correlation"
+    print(corr, title = "\nCorrelation:", ...)
+  }
+
+  cat("\nRandom effects:\n")
+  print(x$confint_ranef, digits = digits)
+
+  cat("\nVariance model:\n")
+  print(x$confint_errmod, digits = digits)
+
+  if (x$transformations == "mkin") {
+    cat("\nBacktransformed parameters:\n")
+    print(x$confint_back, digits = digits)
+  }
+
+  printSFORB <- !is.null(x$SFORB)
+  if(printSFORB){
+    cat("\nEstimated Eigenvalues of SFORB model(s):\n")
+    print(x$SFORB, digits = digits,...)
+  }
+
+  printff <- !is.null(x$ff)
+  if(printff){
+    cat("\nResulting formation fractions:\n")
+    print(data.frame(ff = x$ff), digits = digits,...)
+  }
+
+  printdistimes <- !is.null(x$distimes)
+  if(printdistimes){
+    cat("\nEstimated disappearance times:\n")
+    print(x$distimes, digits = digits,...)
+  }
+
+  if (x$print_data){
+    cat("\nData:\n")
+    print(format(x$data, digits = digits, ...), row.names = FALSE)
+  }
+
+  invisible(x)
+}
-- 
cgit v1.2.1


From 28c0dff7d7191f854be610b5384e965d9b191f98 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 24 Feb 2021 14:46:10 +0100
Subject: Reset graphical parameters with on.exit()

plot.mixed.mmkin did not reset graphical parameters at all. The other
plotting functions did not use on.exit, so this change should make the
use of the plotting functions safer.
---
 R/mkinparplot.R      | 4 ++--
 R/plot.mixed.mmkin.R | 1 +
 R/plot.mkinfit.R     | 2 +-
 R/plot.mmkin.R       | 3 +--
 4 files changed, 5 insertions(+), 5 deletions(-)

(limited to 'R')

diff --git a/R/mkinparplot.R b/R/mkinparplot.R
index f9abab5b..8cae30fb 100644
--- a/R/mkinparplot.R
+++ b/R/mkinparplot.R
@@ -32,7 +32,8 @@ mkinparplot <- function(object) {
               fractions.optim = length(fractions.optim))
   n.plot <- n.plot[n.plot > 0]
 
-  oldpars <- par(no.readonly = TRUE)
+  oldpar <- par(no.readonly = TRUE)
+  on.exit(par(oldpar, no.readonly = TRUE))
   layout(matrix(1:length(n.plot), ncol = 1), heights = n.plot + 1)
 
   s <- summary(object)
@@ -71,5 +72,4 @@ mkinparplot <- function(object) {
 	   as.numeric(values.upper.nonInf), parname_index,
 	   code = 3, angle = 90, length = 0.05))
   }
-  par(oldpars)
 }
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 1674d855..21399496 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -167,6 +167,7 @@ plot.mixed.mmkin <- function(x,
 
   # Start of graphical section
   oldpar <- par(no.readonly = TRUE)
+  on.exit(par(oldpar, no.readonly = TRUE))
 
   n_plot_rows = length(obs_vars)
   n_plots = n_plot_rows * 2
diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
index eced40a4..2e319aae 100644
--- a/R/plot.mkinfit.R
+++ b/R/plot.mkinfit.R
@@ -161,6 +161,7 @@ plot.mkinfit <- function(x, fit = x,
   if (do_layout) {
     # Layout should be restored afterwards
     oldpar <- par(no.readonly = TRUE)
+    on.exit(par(oldpar, no.readonly = TRUE))
 
     # If the observed variables are shown separately, or if requested, do row layout
     if (sep_obs | row_layout) {
@@ -287,7 +288,6 @@ plot.mkinfit <- function(x, fit = x,
         legend = FALSE, frame = frame)
     }
   }
-  if (do_layout) par(oldpar, no.readonly = TRUE)
 }
 
 #' @rdname plot.mkinfit
diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R
index f8ed1f9a..2166b30e 100644
--- a/R/plot.mmkin.R
+++ b/R/plot.mmkin.R
@@ -65,6 +65,7 @@ plot.mmkin <- function(x, main = "auto", legends = 1,
 {
 
   oldpar <- par(no.readonly = TRUE)
+  on.exit(par(oldpar, no.readonly = TRUE))
 
   n.m <- nrow(x)
   n.d <- ncol(x)
@@ -153,6 +154,4 @@ plot.mmkin <- function(x, main = "auto", legends = 1,
     }
     mtext(paste(fit_name, "residuals"), cex = cex, line = 0.4)
   }
-
-  par(oldpar, no.readonly = TRUE)
 }
-- 
cgit v1.2.1


From c73b2f30ec836c949885784ab576e814eb8070a9 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 9 Mar 2021 17:35:47 +0100
Subject: Some improvements for borderline cases

- fit_with_errors for saem()
- test_log_parms for mean_degparms() and saem()
---
 R/nlme.R       | 37 ++++++++++++++++++++++++++++++++++---
 R/nlme.mmkin.R |  2 +-
 R/saem.R       | 31 ++++++++++++++++++++++++++-----
 3 files changed, 61 insertions(+), 9 deletions(-)

(limited to 'R')

diff --git a/R/nlme.R b/R/nlme.R
index 9215aab0..d235a094 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -36,7 +36,7 @@
 #' nlme_f <- nlme_function(f)
 #' # These assignments are necessary for these objects to be
 #' # visible to nlme and augPred when evaluation is done by
-#' # pkgdown to generated the html docs.
+#' # pkgdown to generate the html docs.
 #' assign("nlme_f", nlme_f, globalenv())
 #' assign("grouped_data", grouped_data, globalenv())
 #'
@@ -130,13 +130,44 @@ nlme_function <- function(object) {
 #'   fixed and random effects, in the format required by the start argument of
 #'   nlme for the case of a single grouping variable ds.
 #' @param random Should a list with fixed and random effects be returned?
+#' @param test_log_parms If TRUE, log parameters are only considered in
+#'   the mean calculations if their untransformed counterparts (most likely
+#'   rate constants) pass the t-test for significant difference from zero.
+#' @param conf.level Possibility to adjust the required confidence level
+#'   for parameter that are tested if requested by 'test_log_parms'.
 #' @export
-mean_degparms <- function(object, random = FALSE) {
+mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
+{
   if (nrow(object) > 1) stop("Only row objects allowed")
   parm_mat_trans <- sapply(object, parms, transformed = TRUE)
+
+  if (test_log_parms) {
+      parm_mat_dim <- dim(parm_mat_trans)
+      parm_mat_dimnames <- dimnames(parm_mat_trans)
+
+      log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE)
+      log_parm_names <- gsub("^log_", "", log_parm_trans_names)
+
+      t_test_back_OK <- matrix(
+        sapply(object, function(o) {
+          suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level))
+        }), nrow = length(log_parm_names))
+      rownames(t_test_back_OK) <- log_parm_trans_names
+
+      parm_mat_trans_OK <- parm_mat_trans
+      for (trans_parm in log_parm_trans_names) {
+        parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ],
+          parm_mat_trans[trans_parm, ], NA)
+      }
+    } else {
+    parm_mat_trans_OK <- parm_mat_trans
+  }
+
   mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms))
   degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE]
-  fixed <- apply(degparm_mat_trans, 1, mean)
+  degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE]
+
+  fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE)
   if (random) {
     random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed))
     # If we only have one parameter, apply returns a vector so we get a single row
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index ff1f2fff..306600c6 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -24,7 +24,7 @@ get_deg_func <- function() {
 #' This functions sets up a nonlinear mixed effects model for an mmkin row
 #' object. An mmkin row object is essentially a list of mkinfit objects that
 #' have been obtained by fitting the same model to a list of datasets.
-#' 
+#'
 #' Note that the convergence of the nlme algorithms depends on the quality
 #' of the data. In degradation kinetics, we often only have few datasets
 #' (e.g. data for few soils) and complicated degradation models, which may
diff --git a/R/saem.R b/R/saem.R
index fd2a77b4..460edede 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -24,8 +24,16 @@ utils::globalVariables(c("predicted", "std"))
 #'   SFO or DFOP is used for the parent and there is either no metabolite or one.
 #' @param degparms_start Parameter values given as a named numeric vector will
 #'   be used to override the starting values obtained from the 'mmkin' object.
+#' @param test_log_parms If TRUE, an attempt is made to use more robust starting
+#'   values for population parameters fitted as log parameters in mkin (like
+#'   rate constants) by only considering rate constants that pass the t-test
+#'   when calculating mean degradation parameters using [mean_degparms].
+#' @param conf.level Possibility to adjust the required confidence level
+#'   for parameter that are tested if requested by 'test_log_parms'.
 #' @param solution_type Possibility to specify the solution type in case the
 #'   automatic choice is not desired
+#' @param fail_with_errors Should a failure to compute standard errors
+#'   from the inverse of the Fisher Information Matrix be a failure?
 #' @param quiet Should we suppress the messages saemix prints at the beginning
 #'   and the end of the optimisation process?
 #' @param control Passed to [saemix::saemix]
@@ -51,7 +59,7 @@ utils::globalVariables(c("predicted", "std"))
 #' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
 #' # functions from saemix
 #' library(saemix)
-#' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so))
+#' compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)
 #' plot(f_saem_fomc$so, plot.type = "convergence")
 #' plot(f_saem_fomc$so, plot.type = "individual.fit")
 #' plot(f_saem_fomc$so, plot.type = "npde")
@@ -59,7 +67,7 @@ utils::globalVariables(c("predicted", "std"))
 #'
 #' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
 #' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
-#' compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
+#' compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so)
 #'
 #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
 #'   A1 = mkinsub("SFO"))
@@ -104,19 +112,32 @@ saem <- function(object, ...) UseMethod("saem")
 saem.mmkin <- function(object,
   transformations = c("mkin", "saemix"),
   degparms_start = numeric(),
+  test_log_parms = FALSE,
+  conf.level = 0.6,
   solution_type = "auto",
   control = list(displayProgress = FALSE, print = FALSE,
     save = FALSE, save.graphs = FALSE),
+  fail_with_errors = TRUE,
   verbose = FALSE, quiet = FALSE, ...)
 {
   transformations <- match.arg(transformations)
   m_saemix <- saemix_model(object, verbose = verbose,
-    degparms_start = degparms_start, solution_type = solution_type,
+    degparms_start = degparms_start,
+    test_log_parms = test_log_parms, conf.level = conf.level,
+    solution_type = solution_type,
     transformations = transformations, ...)
   d_saemix <- saemix_data(object, verbose = verbose)
 
   fit_time <- system.time({
     utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
+    FIM_failed <- NULL
+    if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects")
+    if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) {
+      FIM_failed <- c(FIM_failed, "random effects and residual error parameters")
+    }
+    if (!is.null(FIM_failed) & fail_with_errors) {
+      stop("Could not invert FIM for ", paste(FIM_failed, collapse = " and "))
+    }
   })
 
   transparms_optim <- f_saemix@results@fixed.effects
@@ -203,13 +224,13 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
 #' @return An [saemix::SaemixModel] object.
 #' @export
 saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
-  degparms_start = numeric(), verbose = FALSE, ...)
+  degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...)
 {
   if (nrow(object) > 1) stop("Only row objects allowed")
 
   mkin_model <- object[[1]]$mkinmod
 
-  degparms_optim <-  mean_degparms(object)
+  degparms_optim <-  mean_degparms(object, test_log_parms = test_log_parms)
   if (transformations == "saemix") {
     degparms_optim <- backtransform_odeparms(degparms_optim,
       object[[1]]$mkinmod,
-- 
cgit v1.2.1


From cb112e53163f9dc63d439dba50ca051877d67a79 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 16 Mar 2021 16:47:17 +0100
Subject: Convenience option to set nbiter.saemix

---
 R/saem.R | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/saem.R b/R/saem.R
index 460edede..184890f4 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -36,7 +36,9 @@ utils::globalVariables(c("predicted", "std"))
 #'   from the inverse of the Fisher Information Matrix be a failure?
 #' @param quiet Should we suppress the messages saemix prints at the beginning
 #'   and the end of the optimisation process?
-#' @param control Passed to [saemix::saemix]
+#' @param nbiter.saemix Convenience option to increase the number of
+#'   iterations
+#' @param control Passed to [saemix::saemix].
 #' @param \dots Further parameters passed to [saemix::saemixModel].
 #' @return An S3 object of class 'saem.mmkin', containing the fitted
 #'   [saemix::SaemixObject] as a list component named 'so'. The
@@ -115,7 +117,9 @@ saem.mmkin <- function(object,
   test_log_parms = FALSE,
   conf.level = 0.6,
   solution_type = "auto",
+  nbiter.saemix = c(300, 100),
   control = list(displayProgress = FALSE, print = FALSE,
+    nbiter.saemix = nbiter.saemix,
     save = FALSE, save.graphs = FALSE),
   fail_with_errors = TRUE,
   verbose = FALSE, quiet = FALSE, ...)
-- 
cgit v1.2.1


From 6d6dc7d53bf99b088af3488588574afc832fb7fe Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 19 Mar 2021 11:22:07 +0100
Subject: test_log_parms for plot.mixed.mmkin, roxygen run

---
 R/mixed.mmkin.R      | 3 ++-
 R/plot.mixed.mmkin.R | 8 +++++++-
 2 files changed, 9 insertions(+), 2 deletions(-)

(limited to 'R')

diff --git a/R/mixed.mmkin.R b/R/mixed.mmkin.R
index 7aa5edd5..682a9a34 100644
--- a/R/mixed.mmkin.R
+++ b/R/mixed.mmkin.R
@@ -3,6 +3,8 @@
 #' @param object An [mmkin] row object
 #' @param method The method to be used
 #' @param \dots Currently not used
+#' @return An object of class 'mixed.mmkin' which has the observed data in a
+#'   single dataframe which is convenient for plotting
 #' @examples
 #' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
 #' n_biphasic <- 8
@@ -54,7 +56,6 @@ mixed.mmkin <- function(object, method = c("none"), ...) {
   if (nrow(object) > 1) stop("Only row objects allowed")
 
   method <- match.arg(method)
-  if (method == "default") method = c("naive", "nlme")
 
   ds_names <- colnames(object)
   res <- list(mmkin = object, mkinmod = object[[1]]$mkinmod)
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 21399496..f0682c10 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -10,6 +10,10 @@ utils::globalVariables("ds")
 #'   `resplot = "time"`.
 #' @param pred_over Named list of alternative predictions as obtained
 #'   from [mkinpredict] with a compatible [mkinmod].
+#' @param test_log_parms Passed to [mean_degparms] in the case of an
+#'   [mixed.mmkin] object
+#' @param conf.level Passed to [mean_degparms] in the case of an
+#'   [mixed.mmkin] object
 #' @param rel.height.legend The relative height of the legend shown on top
 #' @param rel.height.bottom The relative height of the bottom plot row
 #' @param ymax Vector of maximum y axis values
@@ -58,6 +62,8 @@ plot.mixed.mmkin <- function(x,
   xlim = range(x$data$time),
   resplot = c("predicted", "time"),
   pred_over = NULL,
+  test_log_parms = FALSE,
+  conf.level = 0.6,
   ymax = "auto", maxabs = "auto",
   ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
   nrow.legend = ceiling((length(i) + 1) / ncol.legend),
@@ -76,7 +82,7 @@ plot.mixed.mmkin <- function(x,
   backtransform = TRUE
 
   if (identical(class(x), "mixed.mmkin")) {
-    degparms_pop <- mean_degparms(x$mmkin)
+    degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms, conf.level = conf.level)
 
     degparms_tmp <- parms(x$mmkin, transformed = TRUE)
     degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ]))
-- 
cgit v1.2.1


From 34d1c5f23edfb60548bc5a9dd99c2f3af92acef1 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Sat, 20 Mar 2021 21:26:40 +0100
Subject: Fix mkin calculation of saemix residuals

---
 R/saem.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/saem.R b/R/saem.R
index 184890f4..6f28a47a 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -164,7 +164,7 @@ saem.mmkin <- function(object,
     xidep = return_data[c("time", "name")])
 
   return_data <- transform(return_data,
-    residual = predicted - value,
+    residual = value - predicted,
     std = sigma_twocomp(predicted,
       f_saemix@results@respar[1], f_saemix@results@respar[2]))
   return_data <- transform(return_data,
-- 
cgit v1.2.1


From c6eb6b2bb598002523c3d34d71b0e4a99671ccd6 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 9 Jun 2021 16:53:31 +0200
Subject: Rudimentary support for setting up nlmixr models

- All degradation models are specified as ODE models. This appears to be
fast enough
- Error models are being translated to nlmixr as close to the mkin error
model as possible. When using the 'saem' backend, it appears not to be
possible to use the same error model for more than one observed variable
- No support yet for models with parallel formation of metabolites, where
the ilr transformation is used in mkin per default
- There is a bug in nlmixr which appears to be triggered if the data are
not balanced, see nlmixrdevelopment/nlmixr#530
- There is a print and a plot method, the summary method is not finished
---
 R/mean_degparms.R        |  61 +++++++
 R/nlme.R                 |  55 ------
 R/nlme.mmkin.R           |   2 +-
 R/nlmixr.R               | 467 +++++++++++++++++++++++++++++++++++++++++++++++
 R/plot.mixed.mmkin.R     |  17 ++
 R/saem.R                 |   6 +
 R/summary.nlmixr.mmkin.R | 250 +++++++++++++++++++++++++
 7 files changed, 802 insertions(+), 56 deletions(-)
 create mode 100644 R/mean_degparms.R
 create mode 100644 R/nlmixr.R
 create mode 100644 R/summary.nlmixr.mmkin.R

(limited to 'R')

diff --git a/R/mean_degparms.R b/R/mean_degparms.R
new file mode 100644
index 00000000..ec7f4342
--- /dev/null
+++ b/R/mean_degparms.R
@@ -0,0 +1,61 @@
+#' Calculate mean degradation parameters for an mmkin row object
+#'
+#' @return If random is FALSE (default), a named vector containing mean values
+#'   of the fitted degradation model parameters. If random is TRUE, a list with
+#'   fixed and random effects, in the format required by the start argument of
+#'   nlme for the case of a single grouping variable ds.
+#' @param random Should a list with fixed and random effects be returned?
+#' @param test_log_parms If TRUE, log parameters are only considered in
+#'   the mean calculations if their untransformed counterparts (most likely
+#'   rate constants) pass the t-test for significant difference from zero.
+#' @param conf.level Possibility to adjust the required confidence level
+#'   for parameter that are tested if requested by 'test_log_parms'.
+#' @export
+mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
+{
+  if (nrow(object) > 1) stop("Only row objects allowed")
+  parm_mat_trans <- sapply(object, parms, transformed = TRUE)
+
+  if (test_log_parms) {
+      parm_mat_dim <- dim(parm_mat_trans)
+      parm_mat_dimnames <- dimnames(parm_mat_trans)
+
+      log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE)
+      log_parm_names <- gsub("^log_", "", log_parm_trans_names)
+
+      t_test_back_OK <- matrix(
+        sapply(object, function(o) {
+          suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level))
+        }), nrow = length(log_parm_names))
+      rownames(t_test_back_OK) <- log_parm_trans_names
+
+      parm_mat_trans_OK <- parm_mat_trans
+      for (trans_parm in log_parm_trans_names) {
+        parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ],
+          parm_mat_trans[trans_parm, ], NA)
+      }
+    } else {
+    parm_mat_trans_OK <- parm_mat_trans
+  }
+
+  mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms))
+  degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE]
+  degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE]
+
+  fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE)
+  if (random) {
+    random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed))
+    # If we only have one parameter, apply returns a vector so we get a single row
+    if (nrow(degparm_mat_trans) == 1) random <- t(random)
+    rownames(random) <- levels(nlme_data(object)$ds)
+
+    # For nlmixr we can specify starting values for standard deviations eta, and
+    # we ignore uncertain parameters if test_log_parms is FALSE
+    eta <- apply(degparm_mat_trans_OK, 1, sd, na.rm = TRUE)
+
+    return(list(fixed = fixed, random = list(ds = random), eta = eta))
+  } else {
+    return(fixed)
+  }
+}
+
diff --git a/R/nlme.R b/R/nlme.R
index d235a094..8762f137 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -124,61 +124,6 @@ nlme_function <- function(object) {
   return(model_function)
 }
 
-#' @rdname nlme
-#' @return If random is FALSE (default), a named vector containing mean values
-#'   of the fitted degradation model parameters. If random is TRUE, a list with
-#'   fixed and random effects, in the format required by the start argument of
-#'   nlme for the case of a single grouping variable ds.
-#' @param random Should a list with fixed and random effects be returned?
-#' @param test_log_parms If TRUE, log parameters are only considered in
-#'   the mean calculations if their untransformed counterparts (most likely
-#'   rate constants) pass the t-test for significant difference from zero.
-#' @param conf.level Possibility to adjust the required confidence level
-#'   for parameter that are tested if requested by 'test_log_parms'.
-#' @export
-mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
-{
-  if (nrow(object) > 1) stop("Only row objects allowed")
-  parm_mat_trans <- sapply(object, parms, transformed = TRUE)
-
-  if (test_log_parms) {
-      parm_mat_dim <- dim(parm_mat_trans)
-      parm_mat_dimnames <- dimnames(parm_mat_trans)
-
-      log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE)
-      log_parm_names <- gsub("^log_", "", log_parm_trans_names)
-
-      t_test_back_OK <- matrix(
-        sapply(object, function(o) {
-          suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level))
-        }), nrow = length(log_parm_names))
-      rownames(t_test_back_OK) <- log_parm_trans_names
-
-      parm_mat_trans_OK <- parm_mat_trans
-      for (trans_parm in log_parm_trans_names) {
-        parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ],
-          parm_mat_trans[trans_parm, ], NA)
-      }
-    } else {
-    parm_mat_trans_OK <- parm_mat_trans
-  }
-
-  mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms))
-  degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE]
-  degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE]
-
-  fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE)
-  if (random) {
-    random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed))
-    # If we only have one parameter, apply returns a vector so we get a single row
-    if (nrow(degparm_mat_trans) == 1) random <- t(random)
-    rownames(random) <- levels(nlme_data(object)$ds)
-    return(list(fixed = fixed, random = list(ds = random)))
-  } else {
-    return(fixed)
-  }
-}
-
 #' @rdname nlme
 #' @importFrom purrr map_dfr
 #' @return A \code{\link{groupedData}} object
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index 306600c6..a1aa32e5 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -135,7 +135,7 @@ nlme.mmkin <- function(model, data = "auto",
     function(el) eval(parse(text = paste(el, 1, sep = "~")))),
   random = pdDiag(fixed),
   groups,
-  start = mean_degparms(model, random = TRUE),
+  start = mean_degparms(model, random = TRUE, test_log_parms = TRUE),
   correlation = NULL, weights = NULL,
   subset, method = c("ML", "REML"),
   na.action = na.fail, naPattern,
diff --git a/R/nlmixr.R b/R/nlmixr.R
new file mode 100644
index 00000000..223b23a1
--- /dev/null
+++ b/R/nlmixr.R
@@ -0,0 +1,467 @@
+utils::globalVariables(c("predicted", "std"))
+
+#' Fit nonlinear mixed models using nlmixr
+#'
+#' This function uses [nlmixr::nlmixr()] as a backend for fitting nonlinear mixed
+#' effects models created from [mmkin] row objects using the Stochastic Approximation
+#' Expectation Maximisation algorithm (SAEM).
+#'
+#' An mmkin row object is essentially a list of mkinfit objects that have been
+#' obtained by fitting the same model to a list of datasets using [mkinfit].
+#'
+#' @importFrom nlmixr nlmixr tableControl
+#' @param object An [mmkin] row object containing several fits of the same
+#'   [mkinmod] model to different datasets
+#' @param est Estimation method passed to [nlmixr::nlmixr]
+#' @param degparms_start Parameter values given as a named numeric vector will
+#'   be used to override the starting values obtained from the 'mmkin' object.
+#' @param test_log_parms If TRUE, an attempt is made to use more robust starting
+#'   values for population parameters fitted as log parameters in mkin (like
+#'   rate constants) by only considering rate constants that pass the t-test
+#'   when calculating mean degradation parameters using [mean_degparms].
+#' @param conf.level Possibility to adjust the required confidence level
+#'   for parameter that are tested if requested by 'test_log_parms'.
+#' @param solution_type Possibility to specify the solution type in case the
+#'   automatic choice is not desired
+#' @param control Passed to [nlmixr::nlmixr].
+#' @param \dots Passed to [nlmixr_model]
+#' @return An S3 object of class 'nlmixr.mmkin', containing the fitted
+#'   [nlmixr::nlmixr] object as a list component named 'nm'. The
+#'   object also inherits from 'mixed.mmkin'.
+#' @seealso [summary.nlmixr.mmkin] [plot.mixed.mmkin]
+#' @examples
+#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
+#'  function(x) subset(x$data[c("name", "time", "value")]))
+#' names(ds) <- paste("Dataset", 6:10)
+#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1)
+#' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc",
+#'   cores = 1, quiet = TRUE)
+#'
+#' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem")
+#' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei")
+#'
+#' f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem")
+#' f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei")
+#'
+#' f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem")
+#' f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei")
+#'
+#' f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem")
+#' f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei")
+#'
+#' f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem")
+#' f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei")
+#'
+#' AIC(
+#'   f_nlmixr_sfo_saem$nm, f_nlmixr_sfo_focei$nm,
+#'   f_nlmixr_fomc_saem$nm, f_nlmixr_fomc_focei$nm,
+#'   f_nlmixr_dfop_saem$nm, f_nlmixr_dfop_focei$nm,
+#'   f_nlmixr_hs_saem$nm, f_nlmixr_hs_focei$nm,
+#'   f_nlmixr_fomc_saem_tc$nm, f_nlmixr_fomc_focei_tc$nm)
+#'
+#' AIC(nlme(f_mmkin_parent["FOMC", ]))
+#' AIC(nlme(f_mmkin_parent["HS", ]))
+#'
+#' # nlme is comparable to nlmixr with focei, saem finds a better
+#' # solution, the two-component error model does not improve it
+#' plot(f_nlmixr_fomc_saem)
+#'
+#' \dontrun{
+#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+#'   A1 = mkinsub("SFO"))
+#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+#'   A1 = mkinsub("SFO"))
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+#'   A1 = mkinsub("SFO"))
+#'
+#' f_mmkin_const <- mmkin(list(
+#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#'   ds, quiet = TRUE, error_model = "const")
+#' f_mmkin_obs <- mmkin(list(
+#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#'   ds, quiet = TRUE, error_model = "obs")
+#' f_mmkin_tc <- mmkin(list(
+#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#'   ds, quiet = TRUE, error_model = "tc")
+#'
+#' # A single constant variance is currently only possible with est = 'focei' in nlmixr
+#' f_nlmixr_sfo_sfo_focei_const <- nlmixr(f_mmkin_const["SFO-SFO", ], est = "focei")
+#' f_nlmixr_fomc_sfo_focei_const <- nlmixr(f_mmkin_const["FOMC-SFO", ], est = "focei")
+#' f_nlmixr_dfop_sfo_focei_const <- nlmixr(f_mmkin_const["DFOP-SFO", ], est = "focei")
+#'
+#' # Variance by variable is supported by 'saem' and 'focei'
+#' f_nlmixr_fomc_sfo_saem_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "saem")
+#' f_nlmixr_fomc_sfo_focei_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "focei")
+#' f_nlmixr_dfop_sfo_saem_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "saem")
+#' f_nlmixr_dfop_sfo_focei_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "focei")
+#'
+#' # Identical two-component error for all variables is only possible with
+#' # est = 'focei' in nlmixr
+#' f_nlmixr_fomc_sfo_focei_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei")
+#' f_nlmixr_dfop_sfo_focei_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei")
+#'
+#' # Two-component error by variable is possible with both estimation methods
+#' # Variance by variable is supported by 'saem' and 'focei'
+#' f_nlmixr_fomc_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "saem",
+#'   error_model = "obs_tc")
+#' f_nlmixr_fomc_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei",
+#'   error_model = "obs_tc")
+#' f_nlmixr_dfop_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "saem",
+#'   error_model = "obs_tc")
+#' f_nlmixr_dfop_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei",
+#'   error_model = "obs_tc")
+#'
+#' AIC(
+#'   f_nlmixr_sfo_sfo_focei_const$nm,
+#'   f_nlmixr_fomc_sfo_focei_const$nm,
+#'   f_nlmixr_dfop_sfo_focei_const$nm,
+#'   f_nlmixr_fomc_sfo_saem_obs$nm,
+#'   f_nlmixr_fomc_sfo_focei_obs$nm,
+#'   f_nlmixr_dfop_sfo_saem_obs$nm,
+#'   f_nlmixr_dfop_sfo_focei_obs$nm,
+#'   f_nlmixr_fomc_sfo_focei_tc$nm,
+#'   f_nlmixr_dfop_sfo_focei_tc$nm,
+#'   f_nlmixr_fomc_sfo_saem_obs_tc$nm,
+#'   f_nlmixr_fomc_sfo_focei_obs_tc$nm,
+#'   f_nlmixr_dfop_sfo_saem_obs_tc$nm,
+#'   f_nlmixr_dfop_sfo_focei_obs_tc$nm
+#' )
+#' # Currently, FOMC-SFO with two-component error by variable fitted by focei gives the
+#' # lowest AIC
+#' plot(f_nlmixr_fomc_sfo_focei_obs_tc)
+#' summary(f_nlmixr_fomc_sfo_focei_obs_tc)
+#' }
+#' @export
+nlmixr.mmkin <- function(object, data = NULL,
+  est = NULL, control = list(),
+  table = tableControl(),
+  error_model = object[[1]]$err_mod,
+  test_log_parms = TRUE,
+  conf.level = 0.6,
+  ...,
+  save = NULL,
+  envir = parent.frame()
+)
+{
+  m_nlmixr <- nlmixr_model(object, est = est,
+    error_model = error_model, add_attributes = TRUE,
+    test_log_parms = test_log_parms, conf.level = conf.level)
+  d_nlmixr <- nlmixr_data(object)
+
+  mean_dp_start <- attr(m_nlmixr, "mean_dp_start")
+  mean_ep_start <- attr(m_nlmixr, "mean_ep_start")
+
+  attributes(m_nlmixr) <- NULL
+
+  fit_time <- system.time({
+    f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est)
+  })
+
+  if (is.null(f_nlmixr$CMT)) {
+    nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "IPRED", "IRES", "IWRES")])
+    nlmixr_df$CMT <- as.character(object[[1]]$data$variable[1])
+  } else {
+    nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "CMT", "IPRED", "IRES", "IWRES")])
+  }
+
+  return_data <- nlmixr_df %>%
+    dplyr::transmute(ds = ID, name = CMT, time = TIME, value = DV,
+      predicted = IPRED, residual = IRES,
+      std = IRES/IWRES, standardized = IWRES)
+
+  bparms_optim <- backtransform_odeparms(f_nlmixr$theta,
+    object[[1]]$mkinmod,
+    object[[1]]$transform_rates,
+    object[[1]]$transform_fractions)
+
+  result <- list(
+    mkinmod = object[[1]]$mkinmod,
+    mmkin = object,
+    transform_rates = object[[1]]$transform_rates,
+    transform_fractions = object[[1]]$transform_fractions,
+    nm = f_nlmixr,
+    est = est,
+    time = fit_time,
+    mean_dp_start = mean_dp_start,
+    mean_ep_start = mean_ep_start,
+    bparms.optim = bparms_optim,
+    bparms.fixed = object[[1]]$bparms.fixed,
+    data = return_data,
+    err_mod = error_model,
+    date.fit = date(),
+    nlmixrversion = as.character(utils::packageVersion("nlmixr")),
+    mkinversion = as.character(utils::packageVersion("mkin")),
+    Rversion = paste(R.version$major, R.version$minor, sep=".")
+  )
+
+  class(result) <- c("nlmixr.mmkin", "mixed.mmkin")
+  return(result)
+}
+
+#' @export
+#' @rdname nlmixr.mmkin
+#' @param x An nlmixr.mmkin object to print
+#' @param digits Number of digits to use for printing
+print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+  cat("Kinetic nonlinear mixed-effects model fit by", x$est, "using nlmixr")
+  cat("\nStructural model:\n")
+  diffs <- x$mmkin[[1]]$mkinmod$diffs
+  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs)
+  writeLines(strwrap(nice_diffs, exdent = 11))
+  cat("\nData:\n")
+  cat(nrow(x$data), "observations of",
+    length(unique(x$data$name)), "variable(s) grouped in",
+    length(unique(x$data$ds)), "datasets\n")
+
+  cat("\nLikelihood:\n")
+  print(data.frame(
+      AIC = AIC(x$nm),
+      BIC = BIC(x$nm),
+      logLik = logLik(x$nm),
+      row.names = " "), digits = digits)
+
+  cat("\nFitted parameters:\n")
+  print(x$nm$parFixed, digits = digits)
+
+  invisible(x)
+}
+
+#' @rdname nlmixr.mmkin
+#' @return An function defining a model suitable for fitting with [nlmixr::nlmixr].
+#' @export
+nlmixr_model <- function(object,
+  est = c("saem", "focei"),
+  degparms_start = "auto",
+  test_log_parms = FALSE, conf.level = 0.6,
+  error_model = object[[1]]$err_mod, add_attributes = FALSE)
+{
+  if (nrow(object) > 1) stop("Only row objects allowed")
+  est = match.arg(est)
+
+  mkin_model <- object[[1]]$mkinmod
+  obs_vars <- names(mkin_model$spec)
+
+  if (error_model == object[[1]]$err_mod) {
+
+    if (length(object[[1]]$mkinmod$spec) > 1 & est == "saem") {
+      if (error_model == "const") {
+        message(
+          "Constant variance for more than one variable is not supported for est = 'saem'\n",
+          "Changing the error model to 'obs' (variance by observed variable)")
+        error_model <- "obs"
+      }
+      if (error_model =="tc") {
+        message(
+          "With est = 'saem', a different error model is required for each observed variable",
+          "Changing the error model to 'obs_tc' (Two-component error for each observed variable)")
+        error_model <- "obs_tc"
+      }
+    }
+  }
+
+  degparms_mmkin <- mean_degparms(object,
+    test_log_parms = test_log_parms,
+    conf.level = conf.level, random = TRUE)
+
+  degparms_optim <- degparms_mmkin$fixed
+
+  degparms_optim <- degparms_mmkin$fixed
+
+  if (degparms_start[1] == "auto") {
+    degparms_start <- degparms_optim
+  }
+  degparms_fixed <- object[[1]]$bparms.fixed
+
+  degparms_optim_back_names <- names(backtransform_odeparms(degparms_optim,
+      object[[1]]$mkinmod,
+      object[[1]]$transform_rates,
+      object[[1]]$transform_fractions))
+  names(degparms_optim_back_names) <- names(degparms_optim)
+
+  odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
+  odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE)
+
+  odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names)
+  odeparms_fixed <- degparms_fixed[odeparms_fixed_names]
+
+  odeini_fixed <- degparms_fixed[odeini_fixed_parm_names]
+  names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names)
+
+  # Definition of the model function
+  f <- function(){}
+
+  ini_block <- "ini({"
+
+  # Initial values for all degradation parameters
+  for (parm_name in names(degparms_optim)) {
+    # As initials for state variables are not transformed,
+    # we need to modify the name here as we want to
+    # use the original name in the model block
+    ini_block <- paste0(
+      ini_block,
+      parm_name, " = ",
+      as.character(degparms_start[parm_name]),
+      "\n",
+      "eta.", parm_name, " ~ ",
+      as.character(degparms_mmkin$eta[parm_name]),
+      "\n"
+    )
+  }
+
+  # Error model parameters
+  error_model_mkin <- object[[1]]$err_mod
+
+  errparm_names_mkin <- names(object[[1]]$errparms)
+  errparms_mkin <- sapply(errparm_names_mkin, function(parm_name) {
+    mean(sapply(object, function(x) x$errparms[parm_name]))
+  })
+
+  sigma_tc_mkin <- errparms_ini <- errparms_mkin[1] +
+    mean(unlist(sapply(object, function(x) x$data$observed)), na.rm = TRUE) *
+      errparms_mkin[2]
+
+  if (error_model == "const") {
+    if (error_model_mkin == "tc") {
+      errparms_ini <- sigma_tc_mkin
+    } else {
+      errparms_ini <- mean(errparms_mkin)
+    }
+    names(errparms_ini) <- "sigma"
+  }
+
+  if (error_model == "obs") {
+    errparms_ini <- switch(error_model_mkin,
+      const = rep(errparms_mkin["sigma"], length(obs_vars)),
+      obs = errparms_mkin,
+      tc = sigma_tc_mkin)
+    names(errparms_ini) <- paste0("sigma_", obs_vars)
+  }
+
+  if (error_model == "tc") {
+    if (error_model_mkin != "tc") {
+      stop("Not supported")
+    } else {
+      errparms_ini <- errparms_mkin
+    }
+  }
+
+  if (error_model == "obs_tc") {
+    if (error_model_mkin != "tc") {
+      stop("Not supported")
+    } else {
+      errparms_ini <- rep(errparms_mkin, length(obs_vars))
+      names(errparms_ini) <- paste0(
+        rep(names(errparms_mkin), length(obs_vars)),
+        "_",
+        rep(obs_vars, each = 2))
+    }
+  }
+
+  for (parm_name in names(errparms_ini)) {
+    ini_block <- paste0(
+      ini_block,
+      parm_name, " = ",
+      as.character(errparms_ini[parm_name]),
+      "\n"
+    )
+  }
+
+  ini_block <- paste0(ini_block, "})")
+
+  body(f)[2] <- parse(text = ini_block)
+
+  model_block <- "model({"
+
+  # Population initial values for the ODE state variables
+  for (parm_name in odeini_optim_parm_names) {
+    model_block <- paste0(
+      model_block,
+      parm_name, "_model = ",
+      parm_name, " + eta.", parm_name, "\n",
+      gsub("(.*)_0", "\\1(0)", parm_name), " = ", parm_name, "_model\n")
+  }
+
+  # Population initial values for log rate constants
+  for (parm_name in grep("^log_", names(degparms_optim), value = TRUE)) {
+    model_block <- paste0(
+      model_block,
+      gsub("^log_", "", parm_name), " = ",
+      "exp(", parm_name, " + eta.", parm_name, ")\n")
+  }
+
+  # Population initial values for logit transformed parameters
+  for (parm_name in grep("_qlogis$", names(degparms_optim), value = TRUE)) {
+    model_block <- paste0(
+      model_block,
+      degparms_optim_back_names[parm_name], " = ",
+      "expit(", parm_name, " + eta.", parm_name, ")\n")
+  }
+
+  # Differential equations
+  model_block <- paste0(
+    model_block,
+    paste(
+      gsub("d_(.*) =", "d/dt(\\1) =", mkin_model$diffs),
+      collapse = "\n"),
+    "\n"
+  )
+
+  # Error model
+  if (error_model == "const") {
+    model_block <- paste0(model_block,
+    paste(paste0(obs_vars, " ~ add(sigma)"), collapse = "\n"))
+  }
+  if (error_model == "obs") {
+    model_block <- paste0(model_block,
+    paste(paste0(obs_vars, " ~ add(sigma_", obs_vars, ")"), collapse = "\n"),
+    "\n")
+  }
+  if (error_model == "tc") {
+    model_block <- paste0(model_block,
+    paste(paste0(obs_vars, " ~ add(sigma_low) + prop(rsd_high)"), collapse = "\n"),
+      "\n")
+  }
+  if (error_model == "obs_tc") {
+    model_block <- paste0(model_block,
+    paste(
+      paste0(obs_vars, " ~ add(sigma_low_", obs_vars, ") + ",
+        "prop(rsd_high_", obs_vars, ")"), collapse = "\n"),
+      "\n")
+  }
+
+  model_block <- paste0(model_block, "})")
+
+  body(f)[3] <- parse(text = model_block)
+
+  if (add_attributes) {
+    attr(f, "mean_dp_start") <- degparms_optim
+    attr(f, "mean_ep_start") <- errparms_ini
+  }
+
+  return(f)
+}
+
+#' @rdname nlmixr.mmkin
+#' @return An dataframe suitable for use with [nlmixr::nlmixr]
+#' @export
+nlmixr_data <- function(object, ...) {
+  if (nrow(object) > 1) stop("Only row objects allowed")
+  d <- lapply(object, function(x) x$data)
+  compartment_map <- 1:length(object[[1]]$mkinmod$spec)
+  names(compartment_map) <- names(object[[1]]$mkinmod$spec)
+  ds_names <- colnames(object)
+
+  ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
+  names(ds_list) <- ds_names
+  ds_nlmixr <- purrr::map_dfr(ds_list, function(x) x, .id = "ds")
+  ds_nlmixr$variable <- as.character(ds_nlmixr$variable)
+  ds_nlmixr_renamed <- data.frame(
+    ID = ds_nlmixr$ds,
+    TIME = ds_nlmixr$time,
+    AMT = 0, EVID = 0,
+    CMT = ds_nlmixr$variable,
+    DV = ds_nlmixr$observed,
+    stringsAsFactors = FALSE)
+
+  return(ds_nlmixr_renamed)
+}
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index f0682c10..1ac62b07 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -40,12 +40,17 @@ utils::globalVariables("ds")
 #'
 #' # For this fit we need to increase pnlsMaxiter, and we increase the
 #' # tolerance in order to speed up the fit for this example evaluation
+#' # It still takes 20 seconds to run
 #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
 #' plot(f_nlme)
 #'
 #' f_saem <- saem(f, transformations = "saemix")
 #' plot(f_saem)
 #'
+#' f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs")
+#' f_nlmix <- nlmix(f_obs)
+#' plot(f_nlmix)
+#'
 #' # We can overlay the two variants if we generate predictions
 #' pred_nlme <- mkinpredict(dfop_sfo,
 #'   f_nlme$bparms.optim[-1],
@@ -109,6 +114,18 @@ plot.mixed.mmkin <- function(x,
     names(degparms_pop) <- degparms_i_names
   }
 
+  if (inherits(x, "nlmixr.mmkin")) {
+    eta_i <- random.effects(x$nm)[-1]
+    names(eta_i) <- gsub("^eta.", "", names(eta_i))
+    degparms_i <- eta_i
+    degparms_pop <- x$nm$theta
+    for (parm_name in names(degparms_i)) {
+      degparms_i[parm_name] <- eta_i[parm_name] + degparms_pop[parm_name]
+    }
+    residual_type = ifelse(standardized, "standardized", "residual")
+    residuals <- x$data[[residual_type]]
+  }
+
   degparms_fixed <- fit_1$fixed$value
   names(degparms_fixed) <- rownames(fit_1$fixed)
   degparms_all <- cbind(as.matrix(degparms_i),
diff --git a/R/saem.R b/R/saem.R
index 6f28a47a..5daf4be8 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -13,6 +13,7 @@ utils::globalVariables(c("predicted", "std"))
 #' psi0 of [saemix::saemixModel()] are the mean values of the parameters found
 #' using [mmkin].
 #'
+#' @importFrom utils packageVersion
 #' @param object An [mmkin] row object containing several fits of the same
 #'   [mkinmod] model to different datasets
 #' @param verbose Should we print information about created objects of
@@ -230,6 +231,11 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
 saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
   degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...)
 {
+  if (packageVersion("saemix") < "3.1.9000") {
+    stop("To use the interface to saemix, you need to install a development version\n",
+      "preferably https://github.com/jranke/saemixextension@warp_combined")
+  }
+
   if (nrow(object) > 1) stop("Only row objects allowed")
 
   mkin_model <- object[[1]]$mkinmod
diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R
new file mode 100644
index 00000000..ae8e32cf
--- /dev/null
+++ b/R/summary.nlmixr.mmkin.R
@@ -0,0 +1,250 @@
+#' Summary method for class "nlmixr.mmkin"
+#'
+#' Lists model equations, initial parameter values, optimised parameters
+#' for fixed effects (population), random effects (deviations from the
+#' population mean) and residual error model, as well as the resulting
+#' endpoints such as formation fractions and DT50 values. Optionally
+#' (default is FALSE), the data are listed in full.
+#'
+#' @param object an object of class [nlmix.mmkin]
+#' @param x an object of class [summary.nlmix.mmkin]
+#' @param data logical, indicating whether the full data should be included in
+#'   the summary.
+#' @param verbose Should the summary be verbose?
+#' @param distimes logical, indicating whether DT50 and DT90 values should be
+#'   included.
+#' @param digits Number of digits to use for printing
+#' @param \dots optional arguments passed to methods like \code{print}.
+#' @return The summary function returns a list obtained in the fit, with at
+#'   least the following additional components
+#'   \item{nlmixrversion, mkinversion, Rversion}{The nlmixr, mkin and R versions used}
+#'   \item{date.fit, date.summary}{The dates where the fit and the summary were
+#'     produced}
+#'   \item{diffs}{The differential equations used in the degradation model}
+#'   \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
+#'   \item{data}{The data}
+#'   \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals}
+#'   \item{confint_back}{Backtransformed parameters, with confidence intervals if available}
+#'   \item{confint_errmod}{Error model parameters with confidence intervals}
+#'   \item{ff}{The estimated formation fractions derived from the fitted
+#'      model.}
+#'   \item{distimes}{The DT50 and DT90 values for each observed variable.}
+#'   \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
+#'   The print method is called for its side effect, i.e. printing the summary.
+#' @importFrom stats predict vcov
+#' @author Johannes Ranke for the mkin specific parts
+#'   nlmixr authors for the parts inherited from nlmixr.
+#' @examples
+#' # Generate five datasets following DFOP-SFO kinetics
+#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"),
+#'  m1 = mkinsub("SFO"), quiet = TRUE)
+#' set.seed(1234)
+#' k1_in <- rlnorm(5, log(0.1), 0.3)
+#' k2_in <- rlnorm(5, log(0.02), 0.3)
+#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3))
+#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3))
+#' k_m1_in <- rlnorm(5, log(0.02), 0.3)
+#'
+#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) {
+#'   mkinpredict(dfop_sfo,
+#'     c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1),
+#'     c(parent = 100, m1 = 0),
+#'     sampling_times)
+#' }
+#'
+#' ds_mean_dfop_sfo <- lapply(1:5, function(i) {
+#'   mkinpredict(dfop_sfo,
+#'     c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i],
+#'       f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]),
+#'     c(parent = 100, m1 = 0),
+#'     sampling_times)
+#' })
+#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5)
+#'
+#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) {
+#'   add_err(ds,
+#'     sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
+#'     n = 1)[[1]]
+#' })
+#'
+#' \dontrun{
+#' # Evaluate using mmkin and nlmixr
+#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
+#'   quiet = TRUE, error_model = "tc", cores = 5)
+#' f_saemix_dfop_sfo <- mkin::saem(f_mmkin_dfop_sfo)
+#' f_nlme_dfop_sfo <- mkin::nlme(f_mmkin_dfop_sfo)
+#' f_nlmixr_dfop_sfo_saem <- nlmixr(f_mmkin_dfop_sfo, est = "saem")
+#' # The following takes a very long time but gives
+#' f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei")
+#' AIC(f_nlmixr_dfop_sfo_saem$nm, f_nlmixr_dfop_sfo_focei$nm)
+#' summary(f_nlmixr_dfop_sfo, data = TRUE)
+#' }
+#'
+#' @export
+summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) {
+
+  mod_vars <- names(object$mkinmod$diffs)
+
+  pnames <- names(object$mean_dp_start)
+  np <- length(pnames)
+
+  conf.int <- confint(object$nm)
+  confint_trans <- as.matrix(conf.int[pnames, c(1, 3, 4)])
+  colnames(confint_trans) <- c("est.", "lower", "upper")
+
+  bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
+    object$transform_rates, object$transform_fractions)
+  bpnames <- names(bp)
+
+  # Transform boundaries of CI for one parameter at a time,
+  # with the exception of sets of formation fractions (single fractions are OK).
+  f_names_skip <- character(0)
+  for (box in mod_vars) { # Figure out sets of fractions to skip
+    f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
+    n_paths <- length(f_names)
+    if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
+  }
+
+  confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
+    dimnames = list(bpnames, colnames(confint_trans)))
+  confint_back[, "est."] <- bp
+
+  for (pname in pnames) {
+    if (!pname %in% f_names_skip) {
+      par.lower <- confint_trans[pname, "lower"]
+      par.upper <- confint_trans[pname, "upper"]
+      names(par.lower) <- names(par.upper) <- pname
+      bpl <- backtransform_odeparms(par.lower, object$mkinmod,
+                                            object$transform_rates,
+                                            object$transform_fractions)
+      bpu <- backtransform_odeparms(par.upper, object$mkinmod,
+                                            object$transform_rates,
+                                            object$transform_fractions)
+      confint_back[names(bpl), "lower"] <- bpl
+      confint_back[names(bpu), "upper"] <- bpu
+    }
+  }
+
+  #  Correlation of fixed effects (inspired by summary.nlme)
+  varFix <- vcov(object$nm)
+  stdFix <- sqrt(diag(varFix))
+  object$corFixed <- array(
+    t(varFix/stdFix)/stdFix,
+    dim(varFix),
+    list(pnames, pnames))
+
+  object$confint_back <- confint_back
+
+  object$date.summary = date()
+  object$use_of_ff = object$mkinmod$use_of_ff
+
+  object$diffs <- object$mkinmod$diffs
+  object$print_data <- data # boolean: Should we print the data?
+  predict(object$nm)
+  so_pred <- object$so@results@predictions
+
+  names(object$data)[4] <- "observed" # rename value to observed
+
+  object$verbose <- verbose
+
+  object$fixed <- object$mmkin_orig[[1]]$fixed
+  object$AIC = AIC(object$so)
+  object$BIC = BIC(object$so)
+  object$logLik = logLik(object$so, method = "is")
+
+  ep <- endpoints(object)
+  if (length(ep$ff) != 0)
+    object$ff <- ep$ff
+  if (distimes) object$distimes <- ep$distimes
+  if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB
+  class(object) <- c("summary.saem.mmkin")
+  return(object)
+}
+
+#' @rdname summary.saem.mmkin
+#' @export
+print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
+  cat("saemix version used for fitting:     ", x$saemixversion, "\n")
+  cat("mkin version used for pre-fitting: ", x$mkinversion, "\n")
+  cat("R version used for fitting:        ", x$Rversion, "\n")
+
+  cat("Date of fit:    ", x$date.fit, "\n")
+  cat("Date of summary:", x$date.summary, "\n")
+
+  cat("\nEquations:\n")
+  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
+  writeLines(strwrap(nice_diffs, exdent = 11))
+
+  cat("\nData:\n")
+  cat(nrow(x$data), "observations of",
+    length(unique(x$data$name)), "variable(s) grouped in",
+    length(unique(x$data$ds)), "datasets\n")
+
+  cat("\nModel predictions using solution type", x$solution_type, "\n")
+
+  cat("\nFitted in", x$time[["elapsed"]],  "s using", paste(x$so@options$nbiter.saemix, collapse = ", "), "iterations\n")
+
+  cat("\nVariance model: ")
+  cat(switch(x$err_mod,
+    const = "Constant variance",
+    obs = "Variance unique to each observed variable",
+    tc = "Two-component variance function"), "\n")
+
+  cat("\nMean of starting values for individual parameters:\n")
+  print(x$mean_dp_start, digits = digits)
+
+  cat("\nFixed degradation parameter values:\n")
+  if(length(x$fixed$value) == 0) cat("None\n")
+  else print(x$fixed, digits = digits)
+
+  cat("\nResults:\n\n")
+  cat("Likelihood computed by importance sampling\n")
+  print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
+      row.names = " "), digits = digits)
+
+  cat("\nOptimised parameters:\n")
+  print(x$confint_trans, digits = digits)
+
+  if (nrow(x$confint_trans) > 1) {
+    corr <- x$corFixed
+    class(corr) <- "correlation"
+    print(corr, title = "\nCorrelation:", ...)
+  }
+
+  cat("\nRandom effects:\n")
+  print(x$confint_ranef, digits = digits)
+
+  cat("\nVariance model:\n")
+  print(x$confint_errmod, digits = digits)
+
+  if (x$transformations == "mkin") {
+    cat("\nBacktransformed parameters:\n")
+    print(x$confint_back, digits = digits)
+  }
+
+  printSFORB <- !is.null(x$SFORB)
+  if(printSFORB){
+    cat("\nEstimated Eigenvalues of SFORB model(s):\n")
+    print(x$SFORB, digits = digits,...)
+  }
+
+  printff <- !is.null(x$ff)
+  if(printff){
+    cat("\nResulting formation fractions:\n")
+    print(data.frame(ff = x$ff), digits = digits,...)
+  }
+
+  printdistimes <- !is.null(x$distimes)
+  if(printdistimes){
+    cat("\nEstimated disappearance times:\n")
+    print(x$distimes, digits = digits,...)
+  }
+
+  if (x$print_data){
+    cat("\nData:\n")
+    print(format(x$data, digits = digits, ...), row.names = FALSE)
+  }
+
+  invisible(x)
+}
-- 
cgit v1.2.1


From 0c9b2f0e3c8ce65cb790c9e048476784cbbea070 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 11 Jun 2021 11:14:45 +0200
Subject: Finished 'summary.nlmixr.mmkin', checks, docs

---
 R/endpoints.R            |  4 ++--
 R/mean_degparms.R        |  3 ++-
 R/nlmixr.R               | 29 +++++++++++++++++++++-------
 R/summary.nlmixr.mmkin.R | 50 ++++++++++++++++++++++++------------------------
 4 files changed, 51 insertions(+), 35 deletions(-)

(limited to 'R')

diff --git a/R/endpoints.R b/R/endpoints.R
index f1f47581..6bf52f07 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -10,8 +10,8 @@
 #' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
 #' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
 #'
-#' @param fit An object of class [mkinfit], [nlme.mmkin] or
-#'  [saem.mmkin]. Or another object that has list components
+#' @param fit An object of class [mkinfit], [nlme.mmkin], [saem.mmkin] or
+#'  [nlmixr.mmkin]. Or another object that has list components
 #'  mkinmod containing an [mkinmod] degradation model, and two numeric vectors,
 #'  bparms.optim and bparms.fixed, that contain parameter values
 #'  for that model.
diff --git a/R/mean_degparms.R b/R/mean_degparms.R
index ec7f4342..ec20c068 100644
--- a/R/mean_degparms.R
+++ b/R/mean_degparms.R
@@ -4,6 +4,7 @@
 #'   of the fitted degradation model parameters. If random is TRUE, a list with
 #'   fixed and random effects, in the format required by the start argument of
 #'   nlme for the case of a single grouping variable ds.
+#' @param object An mmkin row object containing several fits of the same model to different datasets
 #' @param random Should a list with fixed and random effects be returned?
 #' @param test_log_parms If TRUE, log parameters are only considered in
 #'   the mean calculations if their untransformed counterparts (most likely
@@ -51,7 +52,7 @@ mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.l
 
     # For nlmixr we can specify starting values for standard deviations eta, and
     # we ignore uncertain parameters if test_log_parms is FALSE
-    eta <- apply(degparm_mat_trans_OK, 1, sd, na.rm = TRUE)
+    eta <- apply(degparm_mat_trans_OK, 1, stats::sd, na.rm = TRUE)
 
     return(list(fixed = fixed, random = list(ds = random), eta = eta))
   } else {
diff --git a/R/nlmixr.R b/R/nlmixr.R
index 223b23a1..98783ca7 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -1,4 +1,7 @@
-utils::globalVariables(c("predicted", "std"))
+utils::globalVariables(c("predicted", "std", "ID", "TIME", "CMT", "DV", "IPRED", "IRES", "IWRES"))
+
+#' @export
+nlmixr::nlmixr
 
 #' Fit nonlinear mixed models using nlmixr
 #'
@@ -10,8 +13,10 @@ utils::globalVariables(c("predicted", "std"))
 #' obtained by fitting the same model to a list of datasets using [mkinfit].
 #'
 #' @importFrom nlmixr nlmixr tableControl
+#' @importFrom dplyr %>%
 #' @param object An [mmkin] row object containing several fits of the same
 #'   [mkinmod] model to different datasets
+#' @param data Not used, the data are extracted from the mmkin row object
 #' @param est Estimation method passed to [nlmixr::nlmixr]
 #' @param degparms_start Parameter values given as a named numeric vector will
 #'   be used to override the starting values obtained from the 'mmkin' object.
@@ -21,22 +26,28 @@ utils::globalVariables(c("predicted", "std"))
 #'   when calculating mean degradation parameters using [mean_degparms].
 #' @param conf.level Possibility to adjust the required confidence level
 #'   for parameter that are tested if requested by 'test_log_parms'.
-#' @param solution_type Possibility to specify the solution type in case the
-#'   automatic choice is not desired
-#' @param control Passed to [nlmixr::nlmixr].
+#' @param data Not used, as the data are extracted from the mmkin row object
+#' @param table Passed to [nlmixr::nlmixr]
+#' @param error_model Possibility to override the error model which is being
+#'   set based on the error model used in the mmkin row object.
+#' @param control Passed to [nlmixr::nlmixr]
 #' @param \dots Passed to [nlmixr_model]
+#' @param save Passed to [nlmixr::nlmixr]
+#' @param envir Passed to [nlmixr::nlmixr]
 #' @return An S3 object of class 'nlmixr.mmkin', containing the fitted
 #'   [nlmixr::nlmixr] object as a list component named 'nm'. The
 #'   object also inherits from 'mixed.mmkin'.
 #' @seealso [summary.nlmixr.mmkin] [plot.mixed.mmkin]
 #' @examples
+#' \dontrun{
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
 #'  function(x) subset(x$data[c("name", "time", "value")]))
 #' names(ds) <- paste("Dataset", 6:10)
+#' 
 #' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1)
 #' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc",
 #'   cores = 1, quiet = TRUE)
-#'
+#' 
 #' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem")
 #' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei")
 #'
@@ -66,7 +77,6 @@ utils::globalVariables(c("predicted", "std"))
 #' # solution, the two-component error model does not improve it
 #' plot(f_nlmixr_fomc_saem)
 #'
-#' \dontrun{
 #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
 #'   A1 = mkinsub("SFO"))
 #' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
@@ -167,7 +177,8 @@ nlmixr.mmkin <- function(object, data = NULL,
   return_data <- nlmixr_df %>%
     dplyr::transmute(ds = ID, name = CMT, time = TIME, value = DV,
       predicted = IPRED, residual = IRES,
-      std = IRES/IWRES, standardized = IWRES)
+      std = IRES/IWRES, standardized = IWRES) %>%
+    dplyr::arrange(ds, name, time)
 
   bparms_optim <- backtransform_odeparms(f_nlmixr$theta,
     object[[1]]$mkinmod,
@@ -227,6 +238,9 @@ print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...)
 }
 
 #' @rdname nlmixr.mmkin
+#' @param add_attributes Should the starting values used for degradation model
+#'   parameters and their distribution and for the error model parameters
+#'   be returned as attributes?
 #' @return An function defining a model suitable for fitting with [nlmixr::nlmixr].
 #' @export
 nlmixr_model <- function(object,
@@ -435,6 +449,7 @@ nlmixr_model <- function(object,
 
   if (add_attributes) {
     attr(f, "mean_dp_start") <- degparms_optim
+    attr(f, "eta_start") <- degparms_mmkin$eta
     attr(f, "mean_ep_start") <- errparms_ini
   }
 
diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R
index ae8e32cf..f2d7c607 100644
--- a/R/summary.nlmixr.mmkin.R
+++ b/R/summary.nlmixr.mmkin.R
@@ -6,8 +6,9 @@
 #' endpoints such as formation fractions and DT50 values. Optionally
 #' (default is FALSE), the data are listed in full.
 #'
-#' @param object an object of class [nlmix.mmkin]
-#' @param x an object of class [summary.nlmix.mmkin]
+#' @importFrom stats confint sd
+#' @param object an object of class [nlmixr.mmkin]
+#' @param x an object of class [summary.nlmixr.mmkin]
 #' @param data logical, indicating whether the full data should be included in
 #'   the summary.
 #' @param verbose Should the summary be verbose?
@@ -23,9 +24,7 @@
 #'   \item{diffs}{The differential equations used in the degradation model}
 #'   \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
 #'   \item{data}{The data}
-#'   \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals}
 #'   \item{confint_back}{Backtransformed parameters, with confidence intervals if available}
-#'   \item{confint_errmod}{Error model parameters with confidence intervals}
 #'   \item{ff}{The estimated formation fractions derived from the fitted
 #'      model.}
 #'   \item{distimes}{The DT50 and DT90 values for each observed variable.}
@@ -78,7 +77,7 @@
 #' # The following takes a very long time but gives
 #' f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei")
 #' AIC(f_nlmixr_dfop_sfo_saem$nm, f_nlmixr_dfop_sfo_focei$nm)
-#' summary(f_nlmixr_dfop_sfo, data = TRUE)
+#' summary(f_nlmixr_dfop_sfo_sfo, data = TRUE)
 #' }
 #'
 #' @export
@@ -134,6 +133,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
     dim(varFix),
     list(pnames, pnames))
 
+  object$confint_trans <- confint_trans
   object$confint_back <- confint_back
 
   object$date.summary = date()
@@ -141,31 +141,29 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
 
   object$diffs <- object$mkinmod$diffs
   object$print_data <- data # boolean: Should we print the data?
-  predict(object$nm)
-  so_pred <- object$so@results@predictions
 
   names(object$data)[4] <- "observed" # rename value to observed
 
   object$verbose <- verbose
 
   object$fixed <- object$mmkin_orig[[1]]$fixed
-  object$AIC = AIC(object$so)
-  object$BIC = BIC(object$so)
-  object$logLik = logLik(object$so, method = "is")
+  object$AIC = AIC(object$nm)
+  object$BIC = BIC(object$nm)
+  object$logLik = logLik(object$nm)
 
   ep <- endpoints(object)
   if (length(ep$ff) != 0)
     object$ff <- ep$ff
   if (distimes) object$distimes <- ep$distimes
   if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB
-  class(object) <- c("summary.saem.mmkin")
+  class(object) <- c("summary.nlmixr.mmkin")
   return(object)
 }
 
-#' @rdname summary.saem.mmkin
+#' @rdname summary.nlmixr.mmkin
 #' @export
-print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
-  cat("saemix version used for fitting:     ", x$saemixversion, "\n")
+print.summary.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
+  cat("nlmixr version used for fitting:   ", x$nlmixrversion, "\n")
   cat("mkin version used for pre-fitting: ", x$mkinversion, "\n")
   cat("R version used for fitting:        ", x$Rversion, "\n")
 
@@ -181,25 +179,29 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
     length(unique(x$data$name)), "variable(s) grouped in",
     length(unique(x$data$ds)), "datasets\n")
 
-  cat("\nModel predictions using solution type", x$solution_type, "\n")
+  cat("\nDegradation model predictions using RxODE\n")
 
-  cat("\nFitted in", x$time[["elapsed"]],  "s using", paste(x$so@options$nbiter.saemix, collapse = ", "), "iterations\n")
+  cat("\nFitted in", x$time[["elapsed"]],  "s\n")
 
   cat("\nVariance model: ")
   cat(switch(x$err_mod,
     const = "Constant variance",
     obs = "Variance unique to each observed variable",
-    tc = "Two-component variance function"), "\n")
+    tc = "Two-component variance function",
+    obs_tc = "Two-component variance unique to each observed variable"), "\n")
 
   cat("\nMean of starting values for individual parameters:\n")
   print(x$mean_dp_start, digits = digits)
 
+  cat("\nMean of starting values for error model parameters:\n")
+  print(x$mean_ep_start, digits = digits)
+
   cat("\nFixed degradation parameter values:\n")
   if(length(x$fixed$value) == 0) cat("None\n")
   else print(x$fixed, digits = digits)
 
   cat("\nResults:\n\n")
-  cat("Likelihood computed by importance sampling\n")
+  cat("Likelihood calculated by", nlmixr::getOfvType(x$nm), " \n")
   print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
       row.names = " "), digits = digits)
 
@@ -212,16 +214,14 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
     print(corr, title = "\nCorrelation:", ...)
   }
 
-  cat("\nRandom effects:\n")
-  print(x$confint_ranef, digits = digits)
+  cat("\nRandom effects (omega):\n")
+  print(x$nm$omega, digits = digits)
 
   cat("\nVariance model:\n")
-  print(x$confint_errmod, digits = digits)
+  print(x$nm$sigma, digits = digits)
 
-  if (x$transformations == "mkin") {
-    cat("\nBacktransformed parameters:\n")
-    print(x$confint_back, digits = digits)
-  }
+  cat("\nBacktransformed parameters:\n")
+  print(x$confint_back, digits = digits)
 
   printSFORB <- !is.null(x$SFORB)
   if(printSFORB){
-- 
cgit v1.2.1


From 88cf130615a6cde0c4e65d14db32fed7f6e43085 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Sat, 12 Jun 2021 11:05:24 +0200
Subject: Small cosmetics

---
 R/dimethenamid_2018.R | 25 +++++++++++++++++++++++++
 R/nlmixr.R            | 20 +++++++++-----------
 2 files changed, 34 insertions(+), 11 deletions(-)

(limited to 'R')

diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 189da618..79018c11 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -18,4 +18,29 @@
 #'   \url{http://registerofquestions.efsa.europa.eu/roqFrontend/outputLoader?output=ON-5211}
 #' @examples
 #' print(dimethenamid_2018)
+#' dmta_ds <- lapply(1:8, function(i) {
+#'   ds_i <- dimethenamid_2018$ds[[i]]$data
+#'   ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
+#'   ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
+#'   ds_i
+#' })
+#' names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+#' dmta_ds[["Borstel"]] <- rbind(dmta_ds[["Borstel 1"]], dmta_ds[["Borstel 2"]])
+#' dmta_ds[["Borstel 1"]] <- NULL
+#' dmta_ds[["Borstel 2"]] <- NULL
+#' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+#' dmta_ds[["Elliot 1"]] <- NULL
+#' dmta_ds[["Elliot 2"]] <- NULL
+#' dfop_sfo3_plus <- mkinmod(
+#'   DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
+#'   M23 = mkinsub("SFO"),
+#'   M27 = mkinsub("SFO"),
+#'   M31 = mkinsub("SFO", "M27", sink = FALSE),
+#'   quiet = TRUE
+#' )
+#' f_dmta_mkin_tc <- mmkin(
+#'   list("DFOP-SFO3+" = dfop_sfo3_plus),
+#'   dmta_ds, quiet = TRUE, error_model = "tc")
+#' nlmixr_model(f_dmta_mkin_tc) # incomplete
+#' # nlmixr(f_dmta_mkin_tc, est = "saem") # not supported (yet)
 "dimethenamid_2018"
diff --git a/R/nlmixr.R b/R/nlmixr.R
index 98783ca7..6e0b5128 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -43,11 +43,11 @@ nlmixr::nlmixr
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
 #'  function(x) subset(x$data[c("name", "time", "value")]))
 #' names(ds) <- paste("Dataset", 6:10)
-#' 
+#'
 #' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1)
 #' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc",
 #'   cores = 1, quiet = TRUE)
-#' 
+#'
 #' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem")
 #' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei")
 #'
@@ -278,20 +278,18 @@ nlmixr_model <- function(object,
     conf.level = conf.level, random = TRUE)
 
   degparms_optim <- degparms_mmkin$fixed
-
-  degparms_optim <- degparms_mmkin$fixed
+  degparms_optim_back <- backtransform_odeparms(degparms_optim,
+      object[[1]]$mkinmod,
+      object[[1]]$transform_rates,
+      object[[1]]$transform_fractions)
+  degparms_optim_back_names <- names(degparms_optim_back)
+  names(degparms_optim_back_names) <- names(degparms_optim)
 
   if (degparms_start[1] == "auto") {
     degparms_start <- degparms_optim
   }
   degparms_fixed <- object[[1]]$bparms.fixed
 
-  degparms_optim_back_names <- names(backtransform_odeparms(degparms_optim,
-      object[[1]]$mkinmod,
-      object[[1]]$transform_rates,
-      object[[1]]$transform_fractions))
-  names(degparms_optim_back_names) <- names(degparms_optim)
-
   odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
   odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE)
 
@@ -307,7 +305,7 @@ nlmixr_model <- function(object,
   ini_block <- "ini({"
 
   # Initial values for all degradation parameters
-  for (parm_name in names(degparms_optim)) {
+  for (parm_name in names(degparms_start)) {
     # As initials for state variables are not transformed,
     # we need to modify the name here as we want to
     # use the original name in the model block
-- 
cgit v1.2.1


From 28197d5fcbaf85b39f4c032b8180d68b6f6a01b3 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 16 Jun 2021 18:03:22 +0200
Subject: Translate formation fractions to nlmixr language

Works for the dimethenamid data, at least for FOCEI. Very little testing
yet. The summary function does not yet handle the new transformations
of formation fractions (that are in fact very old, as they were used
in the very first version of mkin). The test file has no tests yet, just
some code that may be used for testing.
---
 R/dimethenamid_2018.R | 11 +++++--
 R/nlmixr.R            | 79 ++++++++++++++++++++++++++++++++++++++++++++-------
 R/tffm0.R             | 46 ++++++++++++++++++++++++++++++
 3 files changed, 124 insertions(+), 12 deletions(-)
 create mode 100644 R/tffm0.R

(limited to 'R')

diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 79018c11..76b98efe 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -41,6 +41,13 @@
 #' f_dmta_mkin_tc <- mmkin(
 #'   list("DFOP-SFO3+" = dfop_sfo3_plus),
 #'   dmta_ds, quiet = TRUE, error_model = "tc")
-#' nlmixr_model(f_dmta_mkin_tc) # incomplete
-#' # nlmixr(f_dmta_mkin_tc, est = "saem") # not supported (yet)
+#' nlmixr_model(f_dmta_mkin_tc)
+#' f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+#'   control = saemControl(print = 500))
+#' summary(f_dmta_nlmixr_saem)
+#' plot(f_dmta_nlmixr_saem)
+#' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
+#'   control = foceiControl(print = 500))
+#' summary(f_dmta_nlmixr_focei)
+#' plot(f_dmta_nlmixr_focei)
 "dimethenamid_2018"
diff --git a/R/nlmixr.R b/R/nlmixr.R
index 6e0b5128..9c364c4f 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -20,6 +20,9 @@ nlmixr::nlmixr
 #' @param est Estimation method passed to [nlmixr::nlmixr]
 #' @param degparms_start Parameter values given as a named numeric vector will
 #'   be used to override the starting values obtained from the 'mmkin' object.
+#' @param eta_start Standard deviations on the transformed scale given as a
+#'   named numeric vector will be used to override the starting values obtained
+#'   from the 'mmkin' object.
 #' @param test_log_parms If TRUE, an attempt is made to use more robust starting
 #'   values for population parameters fitted as log parameters in mkin (like
 #'   rate constants) by only considering rate constants that pass the t-test
@@ -148,6 +151,8 @@ nlmixr.mmkin <- function(object, data = NULL,
   error_model = object[[1]]$err_mod,
   test_log_parms = TRUE,
   conf.level = 0.6,
+  degparms_start = "auto",
+  eta_start = "auto",
   ...,
   save = NULL,
   envir = parent.frame()
@@ -155,7 +160,9 @@ nlmixr.mmkin <- function(object, data = NULL,
 {
   m_nlmixr <- nlmixr_model(object, est = est,
     error_model = error_model, add_attributes = TRUE,
-    test_log_parms = test_log_parms, conf.level = conf.level)
+    test_log_parms = test_log_parms, conf.level = conf.level,
+    degparms_start = degparms_start, eta_start = eta_start
+  )
   d_nlmixr <- nlmixr_data(object)
 
   mean_dp_start <- attr(m_nlmixr, "mean_dp_start")
@@ -164,7 +171,7 @@ nlmixr.mmkin <- function(object, data = NULL,
   attributes(m_nlmixr) <- NULL
 
   fit_time <- system.time({
-    f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est)
+    f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est, control = control)
   })
 
   if (is.null(f_nlmixr$CMT)) {
@@ -246,7 +253,8 @@ print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...)
 nlmixr_model <- function(object,
   est = c("saem", "focei"),
   degparms_start = "auto",
-  test_log_parms = FALSE, conf.level = 0.6,
+  eta_start = "auto",
+  test_log_parms = TRUE, conf.level = 0.6,
   error_model = object[[1]]$err_mod, add_attributes = FALSE)
 {
   if (nrow(object) > 1) stop("Only row objects allowed")
@@ -278,16 +286,44 @@ nlmixr_model <- function(object,
     conf.level = conf.level, random = TRUE)
 
   degparms_optim <- degparms_mmkin$fixed
+
+  degparms_optim_ilr_names <- grep("^f_.*_ilr", names(degparms_optim), value = TRUE)
+  obs_vars_ilr <- unique(gsub("f_(.*)_ilr.*$", "\\1", degparms_optim_ilr_names))
+  degparms_optim_noilr <- degparms_optim[setdiff(names(degparms_optim),
+    degparms_optim_ilr_names)]
+
   degparms_optim_back <- backtransform_odeparms(degparms_optim,
       object[[1]]$mkinmod,
       object[[1]]$transform_rates,
       object[[1]]$transform_fractions)
-  degparms_optim_back_names <- names(degparms_optim_back)
-  names(degparms_optim_back_names) <- names(degparms_optim)
 
   if (degparms_start[1] == "auto") {
-    degparms_start <- degparms_optim
+    degparms_start <- degparms_optim_noilr
+    for (obs_var_ilr in obs_vars_ilr) {
+      ff_names <- grep(paste0("^f_", obs_var_ilr, "_"),
+        names(degparms_optim_back), value = TRUE)
+      f_tffm0 <- tffm0(degparms_optim_back[ff_names])
+      f_tffm0_qlogis <- qlogis(f_tffm0)
+      names(f_tffm0_qlogis) <- paste0("f_", obs_var_ilr,
+        "_tffm0_", 1:length(f_tffm0), "_qlogis")
+      degparms_start <- c(degparms_start, f_tffm0_qlogis)
+    }
+  }
+
+  if (eta_start[1] == "auto") {
+    eta_start <- degparms_mmkin$eta[setdiff(names(degparms_optim),
+      degparms_optim_ilr_names)]
+    for (obs_var_ilr in obs_vars_ilr) {
+      ff_n <- length(grep(paste0("^f_", obs_var_ilr, "_"),
+        names(degparms_optim_back), value = TRUE))
+      eta_start_ff <- rep(0.3, ff_n)
+      names(eta_start_ff) <- paste0("f_", obs_var_ilr,
+        "_tffm0_", 1:ff_n, "_qlogis")
+      eta_start <- c(eta_start, eta_start_ff)
+    }
   }
+
+
   degparms_fixed <- object[[1]]$bparms.fixed
 
   odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
@@ -315,7 +351,7 @@ nlmixr_model <- function(object,
       as.character(degparms_start[parm_name]),
       "\n",
       "eta.", parm_name, " ~ ",
-      as.character(degparms_mmkin$eta[parm_name]),
+      as.character(eta_start[parm_name]),
       "\n"
     )
   }
@@ -394,7 +430,7 @@ nlmixr_model <- function(object,
   }
 
   # Population initial values for log rate constants
-  for (parm_name in grep("^log_", names(degparms_optim), value = TRUE)) {
+  for (parm_name in grep("^log_", names(degparms_start), value = TRUE)) {
     model_block <- paste0(
       model_block,
       gsub("^log_", "", parm_name), " = ",
@@ -402,13 +438,36 @@ nlmixr_model <- function(object,
   }
 
   # Population initial values for logit transformed parameters
-  for (parm_name in grep("_qlogis$", names(degparms_optim), value = TRUE)) {
+  for (parm_name in grep("_qlogis$", names(degparms_start), value = TRUE)) {
     model_block <- paste0(
       model_block,
-      degparms_optim_back_names[parm_name], " = ",
+      gsub("_qlogis$", "", parm_name), " = ",
       "expit(", parm_name, " + eta.", parm_name, ")\n")
   }
 
+  # Calculate formation fractions from tffm0 transformed values
+  for (obs_var_ilr in obs_vars_ilr) {
+    ff_names <- grep(paste0("^f_", obs_var_ilr, "_"),
+      names(degparms_optim_back), value = TRUE)
+    pattern <- paste0("^f_", obs_var_ilr, "_to_(.*)$")
+    target_vars <- gsub(pattern, "\\1",
+      grep(paste0("^f_", obs_var_ilr, "_to_"), names(degparms_optim_back), value = TRUE))
+    for (i in 1:length(target_vars)) {
+      ff_name <- ff_names[i]
+      ff_line <- paste0(ff_name, " = f_", obs_var_ilr, "_tffm0_", i)
+      if (i > 1) {
+        for (j in (i - 1):1) {
+          ff_line <- paste0(ff_line, " * (1 - f_", obs_var_ilr, "_tffm0_", j , ")")
+        }
+      }
+      model_block <- paste0(
+        model_block,
+        ff_line,
+        "\n"
+      )
+    }
+  }
+
   # Differential equations
   model_block <- paste0(
     model_block,
diff --git a/R/tffm0.R b/R/tffm0.R
new file mode 100644
index 00000000..25787962
--- /dev/null
+++ b/R/tffm0.R
@@ -0,0 +1,46 @@
+#' Transform formation fractions as in the first published mkin version
+#'
+#' The transformed fractions can be restricted between 0 and 1 in model
+#' optimisations. Therefore this transformation was used originally in mkin. It
+#' was later replaced by the [ilr] transformation because the ilr transformed
+#' fractions can assumed to follow normal distribution. As the ilr
+#' transformation is not available in [RxODE] and can therefore not be used in
+#' the nlmixr modelling language, this transformation is currently used for
+#' translating mkin models with formation fractions to more than one target
+#' compartment for fitting with nlmixr in [nlmixr_model]. However,
+#' this implementation cannot be used there, as it is not accessible
+#' from RxODE.
+#'
+#' @param ff Vector of untransformed formation fractions. The sum
+#'   must be smaller or equal to one
+#' @param ff_trans
+#' @return A vector of the transformed formation fractions
+#' @export
+#' @examples
+#' ff_example <- c(
+#'   0.10983681, 0.09035905, 0.08399383
+#' )
+#' ff_example_trans <- tffm0(ff_example)
+#' invtffm0(ff_example_trans)
+tffm0 <- function(ff) {
+  n <- length(ff)
+  res <- numeric(n)
+  f_remaining <- 1
+  for (i in 1:n) {
+    res[i] <- ff[i]/f_remaining
+    f_remaining <- f_remaining - ff[i]
+  }
+  return(res)
+}
+#' @rdname tffm0
+#' @return
+invtffm0 <- function(ff_trans) {
+  n <- length(ff_trans)
+  res <- numeric(n)
+  f_remaining <- 1
+  for (i in 1:n) {
+    res[i] <- ff_trans[i] * f_remaining
+    f_remaining <- f_remaining - res[i]
+  }
+  return(res)
+}
-- 
cgit v1.2.1


From 05baf3bf92cba127fd2319b779db78be86170e5e Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Thu, 17 Jun 2021 13:58:34 +0200
Subject: Let backtransform_odeparms handle nlmixr formation fractions

Also adapt summary.nlmixr.mmkin to correctly handle the way
formation fractions are translated to nlmixr
---
 R/dimethenamid_2018.R    | 14 +++++++++-----
 R/summary.nlmixr.mmkin.R | 14 +++++++-------
 R/tffm0.R                |  6 ++++--
 R/transform_odeparms.R   | 13 +++++++++----
 4 files changed, 29 insertions(+), 18 deletions(-)

(limited to 'R')

diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 76b98efe..6e0bda0c 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -31,6 +31,7 @@
 #' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
 #' dmta_ds[["Elliot 1"]] <- NULL
 #' dmta_ds[["Elliot 2"]] <- NULL
+#' \dontrun{
 #' dfop_sfo3_plus <- mkinmod(
 #'   DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
 #'   M23 = mkinsub("SFO"),
@@ -42,12 +43,15 @@
 #'   list("DFOP-SFO3+" = dfop_sfo3_plus),
 #'   dmta_ds, quiet = TRUE, error_model = "tc")
 #' nlmixr_model(f_dmta_mkin_tc)
-#' f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
-#'   control = saemControl(print = 500))
-#' summary(f_dmta_nlmixr_saem)
-#' plot(f_dmta_nlmixr_saem)
 #' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
-#'   control = foceiControl(print = 500))
+#'   control = nlmixr::foceiControl(print = 500))
 #' summary(f_dmta_nlmixr_focei)
 #' plot(f_dmta_nlmixr_focei)
+#' # saem has a problem with this model/data combination, maybe because of the
+#' # overparameterised error model, to be investigated
+#' #f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+#' #  control = saemControl(print = 500))
+#' #summary(f_dmta_nlmixr_saem)
+#' #plot(f_dmta_nlmixr_saem)
+#' }
 "dimethenamid_2018"
diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R
index f2d7c607..a023f319 100644
--- a/R/summary.nlmixr.mmkin.R
+++ b/R/summary.nlmixr.mmkin.R
@@ -85,11 +85,11 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
 
   mod_vars <- names(object$mkinmod$diffs)
 
-  pnames <- names(object$mean_dp_start)
-  np <- length(pnames)
-
   conf.int <- confint(object$nm)
-  confint_trans <- as.matrix(conf.int[pnames, c(1, 3, 4)])
+  dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start))
+  ndp <- length(dpnames)
+
+  confint_trans <- as.matrix(conf.int[dpnames, c(1, 3, 4)])
   colnames(confint_trans) <- c("est.", "lower", "upper")
 
   bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
@@ -100,7 +100,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
   # with the exception of sets of formation fractions (single fractions are OK).
   f_names_skip <- character(0)
   for (box in mod_vars) { # Figure out sets of fractions to skip
-    f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
+    f_names <- grep(paste("^f", box, sep = "_"), dpnames, value = TRUE)
     n_paths <- length(f_names)
     if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
   }
@@ -109,7 +109,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
     dimnames = list(bpnames, colnames(confint_trans)))
   confint_back[, "est."] <- bp
 
-  for (pname in pnames) {
+  for (pname in dpnames) {
     if (!pname %in% f_names_skip) {
       par.lower <- confint_trans[pname, "lower"]
       par.upper <- confint_trans[pname, "upper"]
@@ -131,7 +131,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes
   object$corFixed <- array(
     t(varFix/stdFix)/stdFix,
     dim(varFix),
-    list(pnames, pnames))
+    list(dpnames, dpnames))
 
   object$confint_trans <- confint_trans
   object$confint_back <- confint_back
diff --git a/R/tffm0.R b/R/tffm0.R
index 25787962..bb5f4cf5 100644
--- a/R/tffm0.R
+++ b/R/tffm0.R
@@ -13,7 +13,8 @@
 #'
 #' @param ff Vector of untransformed formation fractions. The sum
 #'   must be smaller or equal to one
-#' @param ff_trans
+#' @param ff_trans Vector of transformed formation fractions that can be
+#'   restricted to the interval from 0 to 1
 #' @return A vector of the transformed formation fractions
 #' @export
 #' @examples
@@ -33,7 +34,8 @@ tffm0 <- function(ff) {
   return(res)
 }
 #' @rdname tffm0
-#' @return
+#' @export
+#' @return A vector of backtransformed formation fractions for natural use in degradation models
 invtffm0 <- function(ff_trans) {
   n <- length(ff_trans)
   res <- numeric(n)
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 4fe4e5c2..174e7c2d 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -229,13 +229,18 @@ backtransform_odeparms <- function(transparms, mkinmod,
     if (length(trans_f) > 0) {
       if(transform_fractions) {
         if (any(grepl("qlogis", names(trans_f)))) {
-          parms[f_names] <- plogis(trans_f)
+          f_tmp  <- plogis(trans_f)
+          if (any(grepl("_tffm0_.*_qlogis$", names(f_tmp)))) {
+            parms[f_names] <- invtffm0(f_tmp)
+          } else {
+            parms[f_names] <- f_tmp
+          }
         } else {
-          f <- invilr(trans_f)
+          f_tmp <- invilr(trans_f)
           if (spec[[box]]$sink) {
-            parms[f_names] <- f[1:length(f)-1]
+            parms[f_names] <- f_tmp[1:length(f_tmp)-1]
           } else {
-            parms[f_names] <- f
+            parms[f_names] <- f_tmp
           }
         }
       } else {
-- 
cgit v1.2.1


From d9577db290a7fb8944d9a79af59ae90fc00a3eaa Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 23 Jun 2021 17:01:25 +0200
Subject: Fix documentation of default random effects for nlme.mmkin

---
 R/nlme.mmkin.R | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

(limited to 'R')

diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index a1aa32e5..7049a9a1 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -34,10 +34,9 @@ get_deg_func <- function() {
 #' @param data Ignored, data are taken from the mmkin model
 #' @param fixed Ignored, all degradation parameters fitted in the
 #'   mmkin model are used as fixed parameters
-#' @param random If not specified, correlated random effects are set up
-#'   for all optimised degradation model parameters using the log-Cholesky
-#'   parameterization [nlme::pdLogChol] that is also the default of
-#'   the generic [nlme] method.
+#' @param random If not specified, no correlations between random effects are
+#'   set up for the optimised degradation model parameters. This is
+#'   achieved by using the [nlme::pdDiag] method.
 #' @param groups See the documentation of nlme
 #' @param start If not specified, mean values of the fitted degradation
 #'   parameters taken from the mmkin object are used
-- 
cgit v1.2.1


From 8f015900156981ecc2f1f6a1d5a078277ef3f454 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 23 Jun 2021 17:02:20 +0200
Subject: Test log parameters by default when deriving saemix starting
 parameters

---
 R/saem.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/saem.R b/R/saem.R
index 5daf4be8..9db2c04a 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -115,7 +115,7 @@ saem <- function(object, ...) UseMethod("saem")
 saem.mmkin <- function(object,
   transformations = c("mkin", "saemix"),
   degparms_start = numeric(),
-  test_log_parms = FALSE,
+  test_log_parms = TRUE,
   conf.level = 0.6,
   solution_type = "auto",
   nbiter.saemix = c(300, 100),
-- 
cgit v1.2.1


From 40b78bed232798ecbeb72759cdf8d400ea35b31f Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 23 Jul 2021 13:55:34 +0200
Subject: Some example evaluations of dimethenamid data

Evaluations with nlme, saemix and nlmixr are included
---
 R/dimethenamid_2018.R | 33 ++++++++++++++++++++++++---------
 R/mkinsub.R           |  5 -----
 2 files changed, 24 insertions(+), 14 deletions(-)

(limited to 'R')

diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 6e0bda0c..770649e2 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -15,7 +15,7 @@
 #' @source Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria (2018)
 #'   Renewal Assessment Report Dimethenamid-P Volume 3 - B.8 Environmental fate and behaviour
 #'   Rev. 2 - November 2017
-#'   \url{http://registerofquestions.efsa.europa.eu/roqFrontend/outputLoader?output=ON-5211}
+#'   \url{https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716}
 #' @examples
 #' print(dimethenamid_2018)
 #' dmta_ds <- lapply(1:8, function(i) {
@@ -43,15 +43,30 @@
 #'   list("DFOP-SFO3+" = dfop_sfo3_plus),
 #'   dmta_ds, quiet = TRUE, error_model = "tc")
 #' nlmixr_model(f_dmta_mkin_tc)
-#' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
-#'   control = nlmixr::foceiControl(print = 500))
+#' # The focei fit takes about four minutes on my system
+#' system.time(
+#'   f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
+#'     control = nlmixr::foceiControl(print = 500))
+#' )
 #' summary(f_dmta_nlmixr_focei)
 #' plot(f_dmta_nlmixr_focei)
-#' # saem has a problem with this model/data combination, maybe because of the
-#' # overparameterised error model, to be investigated
-#' #f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
-#' #  control = saemControl(print = 500))
-#' #summary(f_dmta_nlmixr_saem)
-#' #plot(f_dmta_nlmixr_saem)
+#' # Using saemix takes about 18 minutes
+#' system.time(
+#'   f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE)
+#' )
+#'
+#' # nlmixr with est = "saem" is pretty fast with default iteration numbers, most
+#' # of the time (about 2.5 minutes) is spent for calculating the log likelihood at the end
+#' # The likelihood calculated for the nlmixr fit is much lower than that found by saemix
+#' # Also, the trace plot and the plot of the individual predictions is not
+#' # convincing for the parent. It seems we are fitting an overparameterised
+#' # model, so the result we get strongly depends on starting parameters and control settings.
+#' system.time(
+#'   f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+#'     control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9))
+#' )
+#' traceplot(f_dmta_nlmixr_saem$nm)
+#' summary(f_dmta_nlmixr_saem)
+#' plot(f_dmta_nlmixr_saem)
 #' }
 "dimethenamid_2018"
diff --git a/R/mkinsub.R b/R/mkinsub.R
index 886f712c..93af3f16 100644
--- a/R/mkinsub.R
+++ b/R/mkinsub.R
@@ -1,8 +1,3 @@
-#' Function to set up a kinetic submodel for one state variable
-#' 
-#' This is a convenience function to set up the lists used as arguments for
-#' \code{\link{mkinmod}}.
-#' 
 #' @rdname mkinmod
 #' @param submodel Character vector of length one to specify the submodel type.
 #'   See \code{\link{mkinmod}} for the list of allowed submodel names.
-- 
cgit v1.2.1


From 5c15ef747568b3a9a9c094b6aa546dc80e3aa87a Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 27 Sep 2021 20:10:01 +0200
Subject: intervals() methods, more DFOP/tc variants

---
 R/intervals.R | 179 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 R/nlmixr.R    |   4 +-
 2 files changed, 181 insertions(+), 2 deletions(-)
 create mode 100644 R/intervals.R

(limited to 'R')

diff --git a/R/intervals.R b/R/intervals.R
new file mode 100644
index 00000000..e2d342f0
--- /dev/null
+++ b/R/intervals.R
@@ -0,0 +1,179 @@
+#' @importFrom nlme intervals
+#' @export
+nlme::intervals
+
+#' Confidence intervals for parameters in saem.mmkin objects
+#'
+#' @param object The fitted saem.mmkin object
+#' @param level The confidence level. Must be the default of 0.95 as this is what
+#'   is available in the saemix object
+#' @param backtransform In case the model was fitted with mkin transformations,
+#'  should we backtransform the parameters where a one to one correlation
+#'  between transformed and backtransformed parameters exists?
+#' @return An object with 'intervals.saem.mmkin' and 'intervals.lme' in the
+#'  class attribute
+#' @export
+intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...)
+{
+  if (!identical(level, 0.95)) {
+    stop("Confidence intervals are only available for a level of 95%")
+  }
+
+  mod_vars <- names(object$mkinmod$diffs)
+
+  pnames <- names(object$mean_dp_start)
+
+  # Confidence intervals are available in the SaemixObject, so
+  # we just need to extract them and put them into a list modelled
+  # after the result of nlme::intervals.lme
+
+  conf.int <- object$so@results@conf.int
+  rownames(conf.int) <- conf.int$name
+  colnames(conf.int)[2] <- "est."
+  confint_trans <- as.matrix(conf.int[pnames, c("lower", "est.", "upper")])
+
+  # Fixed effects
+  # In case objects were produced by earlier versions of saem
+  if (is.null(object$transformations)) object$transformations <- "mkin"
+
+  if (object$transformations == "mkin" & backtransform) {
+    bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
+      object$transform_rates, object$transform_fractions)
+    bpnames <- names(bp)
+
+    # Transform boundaries of CI for one parameter at a time,
+    # with the exception of sets of formation fractions (single fractions are OK).
+    f_names_skip <- character(0)
+    for (box in mod_vars) { # Figure out sets of fractions to skip
+      f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
+      n_paths <- length(f_names)
+      if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
+    }
+
+    confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
+      dimnames = list(bpnames, colnames(confint_trans)))
+    confint_back[, "est."] <- bp
+
+    for (pname in pnames) {
+      if (!pname %in% f_names_skip) {
+        par.lower <- confint_trans[pname, "lower"]
+        par.upper <- confint_trans[pname, "upper"]
+        names(par.lower) <- names(par.upper) <- pname
+        bpl <- backtransform_odeparms(par.lower, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        bpu <- backtransform_odeparms(par.upper, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        confint_back[names(bpl), "lower"] <- bpl
+        confint_back[names(bpu), "upper"] <- bpu
+      }
+    }
+    confint_ret <- confint_back
+  } else {
+    confint_ret <- confint_trans
+  }
+  attr(confint_ret, "label") <- "Fixed effects:"
+
+  # Random effects
+  ranef_ret <- as.matrix(conf.int[paste0("SD.", pnames), c("lower", "est.", "upper")])
+  rownames(ranef_ret) <- paste0(gsub("SD\\.", "sd(", rownames(ranef_ret)), ")")
+  attr(ranef_ret, "label") <- "Random effects:"
+
+
+  # Error model
+  enames <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1")
+  err_ret <- as.matrix(conf.int[enames, c("lower", "est.", "upper")])
+
+  res <- list(
+    fixed = confint_ret,
+    random = ranef_ret,
+    errmod = err_ret
+  )
+  class(res) <- c("intervals.saemix.mmkin", "intervals.lme")
+  attr(res, "level") <- level
+  return(res)
+}
+
+#' Confidence intervals for parameters in nlmixr.mmkin objects
+#'
+#' @param object The fitted saem.mmkin object
+#' @param level The confidence level.
+#' @param backtransform Should we backtransform the parameters where a one to
+#'   one correlation between transformed and backtransformed parameters exists?
+#' @importFrom nlme intervals
+#' @return An object with 'intervals.saem.mmkin' and 'intervals.lme' in the
+#'  class attribute
+#' @export
+intervals.nlmixr.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...)
+{
+
+  # Fixed effects
+  mod_vars <- names(object$mkinmod$diffs)
+
+  conf.int <- confint(object$nm)
+  dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start))
+  ndp <- length(dpnames)
+
+  confint_trans <- as.matrix(conf.int[dpnames, c(3, 1, 4)])
+  colnames(confint_trans) <- c("lower", "est.", "upper")
+
+  if (backtransform) {
+    bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
+      object$transform_rates, object$transform_fractions)
+    bpnames <- names(bp)
+
+    # Transform boundaries of CI for one parameter at a time,
+    # with the exception of sets of formation fractions (single fractions are OK).
+    f_names_skip <- character(0)
+    for (box in mod_vars) { # Figure out sets of fractions to skip
+      f_names <- grep(paste("^f", box, sep = "_"), dpnames, value = TRUE)
+      n_paths <- length(f_names)
+      if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
+    }
+
+    confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
+      dimnames = list(bpnames, colnames(confint_trans)))
+    confint_back[, "est."] <- bp
+
+    for (pname in dpnames) {
+      if (!pname %in% f_names_skip) {
+        par.lower <- confint_trans[pname, "lower"]
+        par.upper <- confint_trans[pname, "upper"]
+        names(par.lower) <- names(par.upper) <- pname
+        bpl <- backtransform_odeparms(par.lower, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        bpu <- backtransform_odeparms(par.upper, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        confint_back[names(bpl), "lower"] <- bpl
+        confint_back[names(bpu), "upper"] <- bpu
+      }
+    }
+    confint_ret <- confint_back
+  } else {
+    confint_ret <- confint_trans
+  }
+  attr(confint_ret, "label") <- "Fixed effects:"
+
+  # Random effects
+  ranef_ret <- as.matrix(data.frame(lower = NA, 
+      "est." = sqrt(diag(object$nm$omega)), upper = NA))
+  rownames(ranef_ret) <- paste0(gsub("eta\\.", "sd(", rownames(ranef_ret)), ")")
+  attr(ranef_ret, "label") <- "Random effects:"
+
+  # Error model
+  enames <- names(object$nm$sigma)
+  err_ret <- as.matrix(conf.int[enames, c(3, 1, 4)])
+  colnames(err_ret) <- c("lower", "est.", "upper")
+
+  res <- list(
+    fixed = confint_ret,
+    random = ranef_ret,
+    errmod = err_ret
+  )
+  class(res) <- c("intervals.nlmixr.mmkin", "intervals.lme")
+  attr(res, "level") <- level
+  return(res)
+}
diff --git a/R/nlmixr.R b/R/nlmixr.R
index 9c364c4f..5f7950ed 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -31,8 +31,8 @@ nlmixr::nlmixr
 #'   for parameter that are tested if requested by 'test_log_parms'.
 #' @param data Not used, as the data are extracted from the mmkin row object
 #' @param table Passed to [nlmixr::nlmixr]
-#' @param error_model Possibility to override the error model which is being
-#'   set based on the error model used in the mmkin row object.
+#' @param error_model Optional argument to override the error model which is
+#'   being set based on the error model used in the mmkin row object.
 #' @param control Passed to [nlmixr::nlmixr]
 #' @param \dots Passed to [nlmixr_model]
 #' @param save Passed to [nlmixr::nlmixr]
-- 
cgit v1.2.1


From cc50f8cad0f608cd2fb9d385f664fc4f53277b2b Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 29 Sep 2021 14:38:38 +0200
Subject: Reduce noise in nlmixr.mmkin output in examples

---
 R/nlmixr.R | 39 +++++++++++++++++++++++++--------------
 1 file changed, 25 insertions(+), 14 deletions(-)

(limited to 'R')

diff --git a/R/nlmixr.R b/R/nlmixr.R
index 5f7950ed..f8fffba9 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -7,7 +7,8 @@ nlmixr::nlmixr
 #'
 #' This function uses [nlmixr::nlmixr()] as a backend for fitting nonlinear mixed
 #' effects models created from [mmkin] row objects using the Stochastic Approximation
-#' Expectation Maximisation algorithm (SAEM).
+#' Expectation Maximisation algorithm (SAEM) or First Order Conditional
+#' Estimation with Interaction (FOCEI).
 #'
 #' An mmkin row object is essentially a list of mkinfit objects that have been
 #' obtained by fitting the same model to a list of datasets using [mkinfit].
@@ -51,20 +52,31 @@ nlmixr::nlmixr
 #' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc",
 #'   cores = 1, quiet = TRUE)
 #'
-#' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem")
-#' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei")
+#' library(nlmixr)
+#' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem",
+#'   control = saemControl(print = 0))
+#' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei",
+#'   control = foceiControl(print = 0))
 #'
-#' f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem")
-#' f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei")
+#' f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem",
+#'   control = saemControl(print = 0))
+#' f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei",
+#'   control = foceiControl(print = 0))
 #'
-#' f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem")
-#' f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei")
+#' f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem",
+#'   control = saemControl(print = 0))
+#' f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei",
+#'   control = foceiControl(print = 0))
 #'
-#' f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem")
-#' f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei")
+#' f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem",
+#'   control = saemControl(print = 0))
+#' f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei",
+#'   control = foceiControl(print = 0))
 #'
-#' f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem")
-#' f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei")
+#' f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem",
+#'   control = saemControl(print = 0))
+#' f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei",
+#'   control = foceiControl(print = 0))
 #'
 #' AIC(
 #'   f_nlmixr_sfo_saem$nm, f_nlmixr_sfo_focei$nm,
@@ -76,9 +88,8 @@ nlmixr::nlmixr
 #' AIC(nlme(f_mmkin_parent["FOMC", ]))
 #' AIC(nlme(f_mmkin_parent["HS", ]))
 #'
-#' # nlme is comparable to nlmixr with focei, saem finds a better
-#' # solution, the two-component error model does not improve it
-#' plot(f_nlmixr_fomc_saem)
+#' # The FOCEI fit of FOMC with constant variance or the tc error model is best
+#' plot(f_nlmixr_fomc_saem_tc)
 #'
 #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
 #'   A1 = mkinsub("SFO"))
-- 
cgit v1.2.1


From 121e38654eea03c91aa117da58fc342feff01756 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Tue, 5 Oct 2021 16:55:41 +0200
Subject: Fix nlmixr for logit transformed parameters

---
 R/nlmixr.R | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

(limited to 'R')

diff --git a/R/nlmixr.R b/R/nlmixr.R
index f8fffba9..55a7bdd7 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -92,11 +92,11 @@ nlmixr::nlmixr
 #' plot(f_nlmixr_fomc_saem_tc)
 #'
 #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
-#'   A1 = mkinsub("SFO"))
+#'   A1 = mkinsub("SFO"), quiet = TRUE)
 #' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
-#'   A1 = mkinsub("SFO"))
+#'   A1 = mkinsub("SFO"), quiet = TRUE)
 #' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
-#'   A1 = mkinsub("SFO"))
+#'   A1 = mkinsub("SFO"), quiet = TRUE)
 #'
 #' f_mmkin_const <- mmkin(list(
 #'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
@@ -108,6 +108,8 @@ nlmixr::nlmixr
 #'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
 #'   ds, quiet = TRUE, error_model = "tc")
 #'
+#' nlmixr_model(f_mmkin_const["SFO-SFO", ])
+#'
 #' # A single constant variance is currently only possible with est = 'focei' in nlmixr
 #' f_nlmixr_sfo_sfo_focei_const <- nlmixr(f_mmkin_const["SFO-SFO", ], est = "focei")
 #' f_nlmixr_fomc_sfo_focei_const <- nlmixr(f_mmkin_const["FOMC-SFO", ], est = "focei")
@@ -450,9 +452,14 @@ nlmixr_model <- function(object,
 
   # Population initial values for logit transformed parameters
   for (parm_name in grep("_qlogis$", names(degparms_start), value = TRUE)) {
+    parm_name_new <- names(
+      backtransform_odeparms(degparms_start[parm_name],
+        object[[1]]$mkinmod,
+        object[[1]]$transform_rates,
+        object[[1]]$transform_fractions))
     model_block <- paste0(
       model_block,
-      gsub("_qlogis$", "", parm_name), " = ",
+      parm_name_new, " = ",
       "expit(", parm_name, " + eta.", parm_name, ")\n")
   }
 
-- 
cgit v1.2.1


From 648637a60a9991cdc8d8d7cbe481cac7f36c953c Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 8 Oct 2021 14:38:19 +0200
Subject: Round initial population values for nlmixr

This avoids numerical instabilities that sometimes occur with the FOCEI
algorithm in nlmixr when the initial values are very close to the
optimum values
---
 R/nlmixr.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

(limited to 'R')

diff --git a/R/nlmixr.R b/R/nlmixr.R
index 55a7bdd7..fd12f555 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -361,10 +361,10 @@ nlmixr_model <- function(object,
     ini_block <- paste0(
       ini_block,
       parm_name, " = ",
-      as.character(degparms_start[parm_name]),
+      as.character(signif(degparms_start[parm_name], 2)),
       "\n",
       "eta.", parm_name, " ~ ",
-      as.character(eta_start[parm_name]),
+      as.character(signif(eta_start[parm_name], 2)),
       "\n"
     )
   }
@@ -422,7 +422,7 @@ nlmixr_model <- function(object,
     ini_block <- paste0(
       ini_block,
       parm_name, " = ",
-      as.character(errparms_ini[parm_name]),
+      as.character(signif(errparms_ini[parm_name], 2)),
       "\n"
     )
   }
-- 
cgit v1.2.1


From d75378911cef79b3ed95daef71bf67db413d2ac8 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Wed, 17 Nov 2021 12:59:49 +0100
Subject: Update required saemix version, update tests

---
 R/saem.R | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

(limited to 'R')

diff --git a/R/saem.R b/R/saem.R
index 9db2c04a..2c20f788 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -231,9 +231,8 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
 saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
   degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...)
 {
-  if (packageVersion("saemix") < "3.1.9000") {
-    stop("To use the interface to saemix, you need to install a development version\n",
-      "preferably https://github.com/jranke/saemixextension@warp_combined")
+  if (packageVersion("saemix") < "3.0") {
+    stop("To use the interface to saemix, you need to install a version >= 3.0\n")
   }
 
   if (nrow(object) > 1) stop("Only row objects allowed")
-- 
cgit v1.2.1