Prepare for v1.0.0v1.0.0

- Improve authorship and copyright information
- Prepare pkgdown config
- Remove dependence on saemix as we need the development version which
is not ready for CRAN
- Temporarily remove saemix interface to check code coverage of the rest

4 files changed, 5 insertions, 806 deletions

diff --git a/R/endpoints.R b/R/endpoints.R
index f1f47581..b5872e68 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] or [nlme.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],
-#'   [summary.nlme.mmkin] and [summary.saem.mmkin].
+#' @note The function is used internally by [summary.mkinfit]
+#'   and [summary.nlme.mmkin]
 #' @author Johannes Ranke
 #' @examples
 #'
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 1674d855..5a0b7412 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], [saem.mmkin] or [nlme.mmkin]
+#' @param x An object of class [mixed.mmkin], [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,15 +39,6 @@ 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,
@@ -91,18 +82,6 @@ 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
deleted file mode 100644
index fd2a77b4..00000000
--- a/R/saem.R
+++ /dev/null
@@ -1,512 +0,0 @@
-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
deleted file mode 100644
index e92c561c..00000000
--- a/R/summary.saem.mmkin.R
+++ /dev/null
@@ -1,268 +0,0 @@
-#' 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)
-}