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

20 files changed, 1209 insertions, 115 deletions

diff --git a/.travis.yml b/.travis.yml
index 6c03b451..60e37230 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -11,6 +11,8 @@ addons:
 cache: packages
 repos:
   CRAN: https://cloud.r-project.org
+r_github_packages:
+  - jranke/saemixextension@warp_combined
 script:
   - R CMD build .
   - R CMD check --no-tests mkin_*.tar.gz
diff --git a/DESCRIPTION b/DESCRIPTION
index 48aaf81f..5b90ef37 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: mkin
 Type: Package
 Title: Kinetic Evaluation of Chemical Degradation Data
-Version: 1.0.4.9000
-Date: 2021-04-21
+Version: 1.0.5
+Date: 2021-06-03
 Authors@R: c(
   person("Johannes", "Ranke", role = c("aut", "cre", "cph"),
     email = "jranke@uni-bremen.de",
@@ -18,11 +18,10 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006,
   note that no warranty is implied for correctness of results or fitness for a
   particular purpose.
 Depends: R (>= 2.15.1), parallel
-Imports: stats, graphics, methods, deSolve, R6, inline (>= 0.3.17), numDeriv,
-  lmtest, pkgbuild, nlme (>= 3.1-151), purrr, saemix (>= 3.1.9000)
+Imports: stats, graphics, methods, deSolve, R6, inline (>= 0.3.19), numDeriv,
+  lmtest, pkgbuild, nlme (>= 3.1-151), purrr, saemix, nlmixr
 Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr,
   benchmarkme, tibble, stats4
-Remotes: github::saemixdevelopment/saemixextension
 License: GPL
 LazyLoad: yes
 LazyData: yes
diff --git a/NAMESPACE b/NAMESPACE
index f2497283..bb4f5f92 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -16,6 +16,7 @@ S3method(mixed,mmkin)
 S3method(mkinpredict,mkinfit)
 S3method(mkinpredict,mkinmod)
 S3method(nlme,mmkin)
+S3method(nlmixr,mmkin)
 S3method(nobs,mkinfit)
 S3method(parms,mkinfit)
 S3method(parms,mmkin)
@@ -30,6 +31,7 @@ S3method(print,mkinmod)
 S3method(print,mmkin)
 S3method(print,nafta)
 S3method(print,nlme.mmkin)
+S3method(print,nlmixr.mmkin)
 S3method(print,saem.mmkin)
 S3method(print,summary.mkinfit)
 S3method(print,summary.nlme.mmkin)
@@ -38,6 +40,7 @@ S3method(residuals,mkinfit)
 S3method(saem,mmkin)
 S3method(summary,mkinfit)
 S3method(summary,nlme.mmkin)
+S3method(summary,nlmixr.mmkin)
 S3method(summary,saem.mmkin)
 S3method(update,mkinfit)
 S3method(update,mmkin)
@@ -86,6 +89,8 @@ export(nafta)
 export(nlme)
 export(nlme_data)
 export(nlme_function)
+export(nlmixr_data)
+export(nlmixr_model)
 export(parms)
 export(plot_err)
 export(plot_res)
@@ -102,6 +107,8 @@ importFrom(R6,R6Class)
 importFrom(grDevices,dev.cur)
 importFrom(lmtest,lrtest)
 importFrom(methods,signature)
+importFrom(nlmixr,nlmixr)
+importFrom(nlmixr,tableControl)
 importFrom(parallel,detectCores)
 importFrom(parallel,mclapply)
 importFrom(parallel,parLapply)
@@ -135,4 +142,5 @@ importFrom(stats,shapiro.test)
 importFrom(stats,update)
 importFrom(stats,vcov)
 importFrom(utils,getFromNamespace)
+importFrom(utils,packageVersion)
 importFrom(utils,write.table)
diff --git a/NEWS.md b/NEWS.md
index d80e152c..03098106 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,18 +1,12 @@
-# mkin 1.0.4.9000
-
-## General
-
-- Switch to a versioning scheme where the fourth version component indicates development versions
+# mkin 1.0.5
 
 ## Mixed-effects models
 
-- Reintroduce the interface to the current development version of saemix, in particular:
-
-- 'saemix_model' and 'saemix_data': Helper functions to set up nonlinear mixed-effects models for mmkin row objects
+- Introduce an interface to nlmixr, supporting estimation methods 'saem' and 'focei': S3 method 'nlmixr.mmkin' using the helper functions 'nlmixr_model' and 'nlmixr_data' to set up nlmixr models for mmkin row objects, with summary and plot methods.
 
-- 'saem': generic function to fit saemix models using 'saemix_model' and 'saemix_data', with a generator 'saem.mmkin', summary and plot methods
+- Reintroduce the interface to current development versions (not on CRAN) of saemix, in particular the generic function 'saem' with a generator 'saem.mmkin', currently using 'saemix_model' and 'saemix_data', summary and plot methods
 
-- 'mean_degparms': New argument 'test_log_parms' that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to obtain more plausible starting parameters for 'saem'
+- 'mean_degparms': New argument 'test_log_parms' that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to obtain more plausible starting parameters for the different mixed-effects model backends
 
 - 'plot.mixed.mmkin': Gains arguments 'test_log_parms' and 'conf.level'
 
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
@@ -125,61 +125,6 @@ nlme_function <- function(object) {
 }
 
 #' @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
 #' @export
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)
+}
diff --git a/build.log b/build.log
index ca1c0481..13f76240 100644
--- a/build.log
+++ b/build.log
@@ -6,5 +6,5 @@
 * creating vignettes ... OK
 * checking for LF line-endings in source and make files and shell scripts
 * checking for empty or unneeded directories
-* building ‘mkin_1.0.4.9000.tar.gz’
+* building ‘mkin_1.0.5.tar.gz’
 
diff --git a/check.log b/check.log
index 7de944a5..f6ee39db 100644
--- a/check.log
+++ b/check.log
@@ -1,19 +1,14 @@
 * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’
-* using R version 4.0.5 (2021-03-31)
+* using R version 4.1.0 (2021-05-18)
 * using platform: x86_64-pc-linux-gnu (64-bit)
 * using session charset: UTF-8
 * using options ‘--no-tests --as-cran’
 * checking for file ‘mkin/DESCRIPTION’ ... OK
 * checking extension type ... Package
-* this is package ‘mkin’ version ‘1.0.4.9000’
+* this is package ‘mkin’ version ‘1.0.5’
 * package encoding: UTF-8
-* checking CRAN incoming feasibility ... NOTE
+* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
 Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
-
-Version contains large components (1.0.4.9000)
-
-Unknown, possibly mis-spelled, fields in DESCRIPTION:
-  ‘Remotes’
 * checking package namespace information ... OK
 * checking package dependencies ... OK
 * checking if this is a source package ... OK
@@ -46,22 +41,68 @@ Unknown, possibly mis-spelled, fields in DESCRIPTION:
 * checking S3 generic/method consistency ... OK
 * checking replacement functions ... OK
 * checking foreign function calls ... OK
-* checking R code for possible problems ... OK
+* checking R code for possible problems ... NOTE
+saemix_model: no visible global function definition for
+  ‘packageVersion’
+Undefined global functions or variables:
+  packageVersion
+Consider adding
+  importFrom("utils", "packageVersion")
+to your NAMESPACE file.
 * checking Rd files ... OK
 * checking Rd metadata ... OK
 * checking Rd line widths ... OK
-* checking Rd cross-references ... OK
+* checking Rd cross-references ... WARNING
+Missing link or links in documentation object 'nlmixr.mmkin.Rd':
+  ‘nlmix_model’ ‘summary.nlmixr.mmkin’
+
+See section 'Cross-references' in the 'Writing R Extensions' manual.
 * checking for missing documentation entries ... OK
 * checking for code/documentation mismatches ... OK
-* checking Rd \usage sections ... OK
+* checking Rd \usage sections ... WARNING
+Undocumented arguments in documentation object 'mean_degparms'
+  ‘object’
+
+Undocumented arguments in documentation object 'nlmixr.mmkin'
+  ‘data’ ‘table’ ‘error_model’ ‘save’ ‘envir’
+Documented arguments not in \usage in documentation object 'nlmixr.mmkin':
+  ‘solution_type’
+
+Functions with \usage entries need to have the appropriate \alias
+entries, and all their arguments documented.
+The \usage entries must correspond to syntactically valid R code.
+See chapter ‘Writing R documentation files’ in the ‘Writing R
+Extensions’ manual.
 * checking Rd contents ... OK
 * checking for unstated dependencies in examples ... OK
 * checking contents of ‘data’ directory ... OK
 * checking data for non-ASCII characters ... OK
+* checking LazyData ... OK
 * checking data for ASCII and uncompressed saves ... OK
 * checking installed files from ‘inst/doc’ ... OK
 * checking files in ‘vignettes’ ... OK
-* checking examples ... OK
+* checking examples ... ERROR
+Running examples in ‘mkin-Ex.R’ failed
+The error most likely occurred in:
+
+> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
+> ### Name: nlmixr.mmkin
+> ### Title: Fit nonlinear mixed models using nlmixr
+> ### Aliases: nlmixr.mmkin print.nlmixr.mmkin nlmixr_model nlmixr_data
+> 
+> ### ** 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")
+Error in nlmixr(f_mmkin_parent["SFO", ], est = "saem") : 
+  could not find function "nlmixr"
+Execution halted
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
@@ -72,9 +113,8 @@ Unknown, possibly mis-spelled, fields in DESCRIPTION:
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: 1 NOTE
+Status: 1 ERROR, 2 WARNINGs, 1 NOTE
 See
   ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
 for details.
 
-
diff --git a/man/mean_degparms.Rd b/man/mean_degparms.Rd
new file mode 100644
index 00000000..92ed4c9d
--- /dev/null
+++ b/man/mean_degparms.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mean_degparms.R
+\name{mean_degparms}
+\alias{mean_degparms}
+\title{Calculate mean degradation parameters for an mmkin row object}
+\usage{
+mean_degparms(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
+}
+\arguments{
+\item{random}{Should a list with fixed and random effects be returned?}
+
+\item{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.}
+
+\item{conf.level}{Possibility to adjust the required confidence level
+for parameter that are tested if requested by 'test_log_parms'.}
+}
+\value{
+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.
+}
+\description{
+Calculate mean degradation parameters for an mmkin row object
+}
diff --git a/man/nlme.Rd b/man/nlme.Rd
index c367868b..e87b7a00 100644
--- a/man/nlme.Rd
+++ b/man/nlme.Rd
@@ -2,36 +2,19 @@
 % Please edit documentation in R/nlme.R
 \name{nlme_function}
 \alias{nlme_function}
-\alias{mean_degparms}
 \alias{nlme_data}
 \title{Helper functions to create nlme models from mmkin row objects}
 \usage{
 nlme_function(object)
 
-mean_degparms(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
-
 nlme_data(object)
 }
 \arguments{
 \item{object}{An mmkin row object containing several fits of the same model to different datasets}
-
-\item{random}{Should a list with fixed and random effects be returned?}
-
-\item{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.}
-
-\item{conf.level}{Possibility to adjust the required confidence level
-for parameter that are tested if requested by 'test_log_parms'.}
 }
 \value{
 A function that can be used with nlme
 
-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.
-
 A \code{\link{groupedData}} object
 }
 \description{
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
index 2fb0488a..a2b45efa 100644
--- a/man/nlme.mmkin.Rd
+++ b/man/nlme.mmkin.Rd
@@ -13,7 +13,7 @@
     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,
diff --git a/man/nlmixr.mmkin.Rd b/man/nlmixr.mmkin.Rd
new file mode 100644
index 00000000..86bbdc9f
--- /dev/null
+++ b/man/nlmixr.mmkin.Rd
@@ -0,0 +1,188 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nlmixr.R
+\name{nlmixr.mmkin}
+\alias{nlmixr.mmkin}
+\alias{print.nlmixr.mmkin}
+\alias{nlmixr_model}
+\alias{nlmixr_data}
+\title{Fit nonlinear mixed models using nlmixr}
+\usage{
+\method{nlmixr}{mmkin}(
+  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()
+)
+
+\method{print}{nlmixr.mmkin}(x, digits = max(3, getOption("digits") - 3), ...)
+
+nlmixr_model(
+  object,
+  est = c("saem", "focei"),
+  degparms_start = "auto",
+  test_log_parms = FALSE,
+  conf.level = 0.6,
+  error_model = object[[1]]$err_mod
+)
+
+nlmixr_data(object, ...)
+}
+\arguments{
+\item{object}{An \link{mmkin} row object containing several fits of the same
+\link{mkinmod} model to different datasets}
+
+\item{est}{Estimation method passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
+
+\item{control}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}.}
+
+\item{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 \link{mean_degparms}.}
+
+\item{conf.level}{Possibility to adjust the required confidence level
+for parameter that are tested if requested by 'test_log_parms'.}
+
+\item{\dots}{Passed to \link{nlmixr_model}}
+
+\item{x}{An nlmixr.mmkin object to print}
+
+\item{digits}{Number of digits to use for printing}
+
+\item{degparms_start}{Parameter values given as a named numeric vector will
+be used to override the starting values obtained from the 'mmkin' object.}
+
+\item{solution_type}{Possibility to specify the solution type in case the
+automatic choice is not desired}
+}
+\value{
+An S3 object of class 'nlmixr.mmkin', containing the fitted
+\link[nlmixr:nlmixr]{nlmixr::nlmixr} object as a list component named 'nm'. The
+object also inherits from 'mixed.mmkin'.
+
+An function defining a model suitable for fitting with \link[nlmixr:nlmixr]{nlmixr::nlmixr}.
+
+An dataframe suitable for use with \link[nlmixr:nlmixr]{nlmixr::nlmixr}
+}
+\description{
+This function uses \code{\link[nlmixr:nlmixr]{nlmixr::nlmixr()}} as a backend for fitting nonlinear mixed
+effects models created from \link{mmkin} row objects using the Stochastic Approximation
+Expectation Maximisation algorithm (SAEM).
+}
+\details{
+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 \link{mkinfit}.
+}
+\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)
+}
+}
+\seealso{
+\link{summary.nlmixr.mmkin} \link{plot.mixed.mmkin}
+}
diff --git a/man/plot.mixed.mmkin.Rd b/man/plot.mixed.mmkin.Rd
index bcab3e74..d87ca22c 100644
--- a/man/plot.mixed.mmkin.Rd
+++ b/man/plot.mixed.mmkin.Rd
@@ -99,12 +99,17 @@ plot(f[, 3:4], standardized = TRUE)
 
 # 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],
diff --git a/man/summary.nlmixr.mmkin.Rd b/man/summary.nlmixr.mmkin.Rd
new file mode 100644
index 00000000..03f0ffb2
--- /dev/null
+++ b/man/summary.nlmixr.mmkin.Rd
@@ -0,0 +1,100 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/summary.nlmixr.mmkin.R
+\name{summary.nlmixr.mmkin}
+\alias{summary.nlmixr.mmkin}
+\title{Summary method for class "nlmixr.mmkin"}
+\usage{
+\method{summary}{nlmixr.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...)
+}
+\arguments{
+\item{object}{an object of class \link{nlmix.mmkin}}
+
+\item{data}{logical, indicating whether the full data should be included in
+the summary.}
+
+\item{verbose}{Should the summary be verbose?}
+
+\item{distimes}{logical, indicating whether DT50 and DT90 values should be
+included.}
+
+\item{\dots}{optional arguments passed to methods like \code{print}.}
+
+\item{x}{an object of class \link{summary.nlmix.mmkin}}
+
+\item{digits}{Number of digits to use for printing}
+}
+\value{
+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.
+}
+\description{
+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.
+}
+\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 = "obs", 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")
+#f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei")
+summary(f_nlmixr_dfop_sfo, data = TRUE)
+}
+
+}
+\author{
+Johannes Ranke for the mkin specific parts
+nlmixr authors for the parts inherited from nlmixr.
+}
diff --git a/man/summary.saem.mmkin.Rd b/man/summary.saem.mmkin.Rd
index 67cb3cbb..86938d31 100644
--- a/man/summary.saem.mmkin.Rd
+++ b/man/summary.saem.mmkin.Rd
@@ -1,30 +1,32 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/summary.saem.mmkin.R
-\name{summary.saem.mmkin}
-\alias{summary.saem.mmkin}
+% Please edit documentation in R/summary.nlmixr.mmkin.R, R/summary.saem.mmkin.R
+\name{print.summary.saem.mmkin}
 \alias{print.summary.saem.mmkin}
+\alias{summary.saem.mmkin}
 \title{Summary method for class "saem.mmkin"}
 \usage{
+\method{print}{summary.saem.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...)
+
 \method{summary}{saem.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...)
 
 \method{print}{summary.saem.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...)
 }
 \arguments{
-\item{object}{an object of class \link{saem.mmkin}}
+\item{x}{an object of class \link{summary.saem.mmkin}}
 
-\item{data}{logical, indicating whether the full data should be included in
-the summary.}
+\item{digits}{Number of digits to use for printing}
 
 \item{verbose}{Should the summary be verbose?}
 
-\item{distimes}{logical, indicating whether DT50 and DT90 values should be
-included.}
-
 \item{\dots}{optional arguments passed to methods like \code{print}.}
 
-\item{x}{an object of class \link{summary.saem.mmkin}}
+\item{object}{an object of class \link{saem.mmkin}}
 
-\item{digits}{Number of digits to use for printing}
+\item{data}{logical, indicating whether the full data should be included in
+the summary.}
+
+\item{distimes}{logical, indicating whether DT50 and DT90 values should be
+included.}
 }
 \value{
 The summary function returns a list based on the \link[saemix:SaemixObject-class]{saemix::SaemixObject}