Improve saem method, add summary

Also make the endpoints function work for saem objects.

13 files changed, 459 insertions, 82 deletions

diff --git a/NAMESPACE b/NAMESPACE
index 3360f861..ed1a46de 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -27,10 +27,12 @@ S3method(print,nafta)
 S3method(print,nlme.mmkin)
 S3method(print,summary.mkinfit)
 S3method(print,summary.nlme.mmkin)
+S3method(print,summary.saem.mmkin)
 S3method(residuals,mkinfit)
 S3method(saem,mmkin)
 S3method(summary,mkinfit)
 S3method(summary,nlme.mmkin)
+S3method(summary,saem.mmkin)
 S3method(update,mkinfit)
 S3method(update,mmkin)
 S3method(update,nlme.mmkin)
@@ -123,5 +125,6 @@ importFrom(stats,residuals)
 importFrom(stats,rnorm)
 importFrom(stats,shapiro.test)
 importFrom(stats,update)
+importFrom(stats,vcov)
 importFrom(utils,getFromNamespace)
 importFrom(utils,write.table)
diff --git a/NEWS.md b/NEWS.md
index 60276151..de421202 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # mkin 0.9.50.4 (unreleased)
 
+- 'saem' generic function to fit saemix models, with a method 'saem.mmkin' and a corresponding 'summary.saem.mmkin'
+
 - 'transform_odeparms', 'backtransform_odeparms': Use logit transformation for solitary fractions like the g parameter of the DFOP model, or formation fractions for a pathway to only one target variable
 
 - 'update' method for 'mmkin' objects
diff --git a/R/endpoints.R b/R/endpoints.R
index e4813db9..f1f47581 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -7,16 +7,22 @@
 #' are equivalent to the rate constants of the DFOP model, but with the
 #' advantage that the SFORB model can also be used for metabolites.
 #'
-#' @param fit An object of class \code{\link{mkinfit}} or
-#'   \code{\link{nlme.mmkin}}
+#' 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
+#'  mkinmod containing an [mkinmod] degradation model, and two numeric vectors,
+#'  bparms.optim and bparms.fixed, that contain parameter values
+#'  for that model.
 #' @importFrom stats optimize
 #' @return A list with a matrix of dissipation times named distimes,
 #'   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 \code{\link{summary.mkinfit}}.
+#' @note The function is used internally by [summary.mkinfit],
+#'   [summary.nlme.mmkin] and [summary.saem.mmkin].
 #' @author Johannes Ranke
-#' @keywords manip
 #' @examples
 #'
 #'   fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
@@ -30,26 +36,9 @@
 #'
 #' @export
 endpoints <- function(fit) {
-  # Calculate dissipation times DT50 and DT90 and formation
-  # fractions as well as SFORB eigenvalues from optimised parameters
-  # 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
   ep <- list()
-  if (inherits(fit, "mkinfit")) {
-    mkinmod <- fit$mkinmod
-    parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
-  } else {
-    if (inherits(fit, "nlme.mmkin")) {
-      mkinmod <- fit$mmkin_orig[[1]]$mkinmod
-      bparms.optim <- backtransform_odeparms(fit$coefficients$fixed,
-        mkinmod,
-        transform_rates = fit$mmkin_orig[[1]]$transform_rates,
-        transform_fractions = fit$mmkin_orig[[1]]$transform_fractions)
-      parms.all <- c(bparms.optim, fit$bparms.fixed)
-    } else {
-      stop("Only implemented for mkinfit and nlme.mmkin objects")
-    }
-  }
+  mkinmod <- fit$mkinmod
+  degparms <- c(fit$bparms.optim, fit$bparms.fixed)
   obs_vars <- names(mkinmod$spec)
   ep$ff <- vector()
   ep$SFORB <- vector()
@@ -61,9 +50,9 @@ endpoints <- function(fit) {
     type = names(mkinmod$map[[obs_var]])[1]
 
     # Get formation fractions if directly fitted, and calculate remaining fraction to sink
-    f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE)
+    f_names = grep(paste("^f", obs_var, sep = "_"), names(degparms), value=TRUE)
     if (length(f_names) > 0) {
-      f_values = parms.all[f_names]
+      f_values = degparms[f_names]
       f_to_sink = 1 - sum(f_values)
       names(f_to_sink) = ifelse(type == "SFORB",
         paste(obs_var, "free", "sink", sep = "_"),
@@ -76,34 +65,34 @@ endpoints <- function(fit) {
 
     # Get the rest
     if (type == "SFO") {
-      k_names = grep(paste("^k", obs_var, sep="_"), names(parms.all), value=TRUE)
-      k_tot = sum(parms.all[k_names])
+      k_names = grep(paste("^k", obs_var, sep="_"), names(degparms), value=TRUE)
+      k_tot = sum(degparms[k_names])
       DT50 = log(2)/k_tot
       DT90 = log(10)/k_tot
       if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) {
         for (k_name in k_names)
         {
-          ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
+          ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
         }
       }
     }
     if (type == "FOMC") {
-      alpha = parms.all["alpha"]
-      beta = parms.all["beta"]
+      alpha = degparms["alpha"]
+      beta = degparms["beta"]
       DT50 = beta * (2^(1/alpha) - 1)
       DT90 = beta * (10^(1/alpha) - 1)
       DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
       ep$distimes[obs_var, c("DT50back")] = DT50_back
     }
     if (type == "IORE") {
-      k_names = grep(paste("^k__iore", obs_var, sep="_"), names(parms.all), value=TRUE)
-      k_tot = sum(parms.all[k_names])
+      k_names = grep(paste("^k__iore", obs_var, sep="_"), names(degparms), value=TRUE)
+      k_tot = sum(degparms[k_names])
       # From the NAFTA kinetics guidance, p. 5
-      n = parms.all[paste("N", obs_var, sep = "_")]
+      n = degparms[paste("N", obs_var, sep = "_")]
       k = k_tot
       # Use the initial concentration of the parent compound
       source_name = mkinmod$map[[1]][[1]]
-      c0 = parms.all[paste(source_name, "0", sep = "_")]
+      c0 = degparms[paste(source_name, "0", sep = "_")]
       alpha = 1 / (n - 1)
       beta = (c0^(1 - n))/(k * (n - 1))
       DT50 = beta * (2^(1/alpha) - 1)
@@ -113,14 +102,14 @@ endpoints <- function(fit) {
       if (mkinmod$use_of_ff == "min") {
         for (k_name in k_names)
         {
-          ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
+          ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
         }
       }
     }
     if (type == "DFOP") {
-      k1 = parms.all["k1"]
-      k2 = parms.all["k2"]
-      g = parms.all["g"]
+      k1 = degparms["k1"]
+      k2 = degparms["k2"]
+      g = degparms["g"]
       f <- function(log_t, x) {
         t <- exp(log_t)
         fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
@@ -144,9 +133,9 @@ endpoints <- function(fit) {
       ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
     }
     if (type == "HS") {
-      k1 = parms.all["k1"]
-      k2 = parms.all["k2"]
-      tb = parms.all["tb"]
+      k1 = degparms["k1"]
+      k2 = degparms["k2"]
+      tb = degparms["tb"]
       DTx <- function(x) {
         DTx.a <- (log(100/(100 - x)))/k1
         DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
@@ -165,11 +154,11 @@ endpoints <- function(fit) {
     }
     if (type == "SFORB") {
       # FOCUS kinetics (2006), p. 60 f
-      k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
+      k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(degparms), value=TRUE)
       k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
-      k_1output = sum(parms.all[k_out_names])
-      k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]
-      k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]
+      k_1output = sum(degparms[k_out_names])
+      k_12 = degparms[paste("k", obs_var, "free", "bound", sep="_")]
+      k_21 = degparms[paste("k", obs_var, "bound", "free", sep="_")]
 
       sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
       b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
@@ -201,7 +190,7 @@ endpoints <- function(fit) {
 
       for (k_out_name in k_out_names)
       {
-        ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output
+        ep$ff[[sub("k_", "", k_out_name)]] = degparms[[k_out_name]] / k_1output
       }
 
       # Return the eigenvalues for comparison with DFOP rate constants
@@ -214,9 +203,9 @@ endpoints <- function(fit) {
     }
     if (type == "logistic") {
       # FOCUS kinetics (2014) p. 67
-      kmax = parms.all["kmax"]
-      k0 = parms.all["k0"]
-      r = parms.all["r"]
+      kmax = degparms["kmax"]
+      k0 = degparms["k0"]
+      r = degparms["r"]
       DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))
       DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))
 
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index 526cb10b..af92e8a1 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -181,14 +181,21 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
   fit_time <- system.time(val <- do.call("nlme.formula", thisCall))
   val$time <- fit_time
 
-  val$mean_dp_start <- mean_dp_start
-  val$mmkin_orig <- model
-  val$data <- thisCall[["data"]]
   val$mkinmod <- model[[1]]$mkinmod
-  val$err_mode <- error_model
+  val$data <- thisCall[["data"]]
+  val$mmkin_orig <- model
+  val$mean_dp_start <- mean_dp_start
   val$transform_rates <- model[[1]]$transform_rates
   val$transform_fractions <- model[[1]]$transform_fractions
   val$solution_type <- model[[1]]$solution_type
+  val$err_mode <- error_model
+
+  val$bparms.optim <- backtransform_odeparms(val$coefficients$fixed,
+    val$mkinmod,
+    transform_rates = val$transform_rates,
+    transform_fractions = val$transform_fractions)
+
+  val$bparms.fixed <- model[[1]]$bparms.fixed
   val$date.fit <- date()
   val$nlmeversion <- as.character(utils::packageVersion("nlme"))
   val$mkinversion <- as.character(utils::packageVersion("mkin"))
diff --git a/R/saemix.R b/R/saemix.R
index 1db8b011..090f0017 100644
--- a/R/saemix.R
+++ b/R/saemix.R
@@ -1,8 +1,8 @@
 #' 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
-#' to the expectation maximisation algorithm (SAEM).
+#' 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].
@@ -23,7 +23,9 @@
 #' @param control Passed to [saemix::saemix]
 #' @param \dots Further parameters passed to [saemix::saemixData]
 #'   and [saemix::saemixModel].
-#' @return An [saemix::SaemixObject].
+#' @return An S3 object of class 'saem.mmkin', containing the fitted
+#'   [saemix::SaemixObject] as a list component named 'so'.
+#' @seealso [summary.saem.mmkin]
 #' @examples
 #' \dontrun{
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
@@ -57,7 +59,7 @@
 #' # Using a single core, the following takes about 6 minutes, using 10 cores
 #' # it is slower instead of faster
 #' f_saem_des <- saem(f_mmkin_des, cores = 1)
-#' compare.saemix(list(f_saemix$so, f_saemix_des$so))
+#' compare.saemix(list(f_saem$so, f_saem_des$so))
 #' }
 #' @export
 saem <- function(object, control, ...) UseMethod("saem")
@@ -79,18 +81,40 @@ saem.mmkin <- function(object,
     tmp <- tempfile()
     grDevices::png(tmp)
   }
-  f_saemix <- saemix::saemix(m_saemix, d_saemix, control)
+  fit_time <- system.time({
+    f_saemix <- saemix::saemix(m_saemix, d_saemix, control)
+    f_saemix <- saemix::saemix.predict(f_saemix)
+  })
   if (suppressPlot) {
     grDevices::dev.off()
     unlink(tmp)
   }
+  transparms_optim = f_saemix@results@fixed.effects
+  names(transparms_optim) = f_saemix@results@name.fixed
+  bparms_optim <- backtransform_odeparms(transparms_optim,
+    object[[1]]$mkinmod,
+    object[[1]]$transform_rates,
+    object[[1]]$transform_fractions)
+
   result <- list(
     mkinmod = object[[1]]$mkinmod,
     mmkin = object,
     solution_type = object[[1]]$solution_type,
     transform_rates = object[[1]]$transform_rates,
     transform_fractions = object[[1]]$transform_fractions,
-    so = f_saemix)
+    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 = nlme_data(object),
+    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) <- "saem.mmkin"
   return(result)
 }
@@ -256,6 +280,7 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
     error.init = error.init,
     verbose = verbose
   )
+  attr(res, "mean_dp_start") <- degparms_optim
   return(res)
 }
 
diff --git a/R/summary.nlme.mmkin.R b/R/summary.nlme.mmkin.R
index 9fdd3f73..7e404e00 100644
--- a/R/summary.nlme.mmkin.R
+++ b/R/summary.nlme.mmkin.R
@@ -75,7 +75,7 @@ summary.nlme.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
   confint_trans <- intervals(object, which = "fixed", level = 1 - alpha)$fixed
   attr(confint_trans, "label") <- NULL
   pnames <- rownames(confint_trans)
-  confint_trans[, "est."]
+
   bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
     object$transform_rates, object$transform_fractions)
   bpnames <- names(bp)
@@ -127,14 +127,12 @@ summary.nlme.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
 
   object$diffs <- object$mkinmod$diffs
   object$print_data <- data
-  if (data) {
-    object$data[["observed"]] <- object$data[["value"]]
-    object$data[["value"]] <- NULL
-    object$data[["predicted"]] <- predict(object)
-    object$data[["residual"]] <- residuals(object, type = "response")
-    object$data[["std"]] <- object$sigma <- 1/attr(object$modelStruct$varStruct, "weights")
-    object$data[["standardized"]] <- residuals(object, type = "pearson")
-  }
+  object$data[["observed"]] <- object$data[["value"]]
+  object$data[["value"]] <- NULL
+  object$data[["predicted"]] <- predict(object)
+  object$data[["residual"]] <- residuals(object, type = "response")
+  object$data[["std"]] <- object$sigma <- 1/attr(object$modelStruct$varStruct, "weights")
+  object$data[["standardized"]] <- residuals(object, type = "pearson")
   object$verbose <- verbose
 
   object$fixed <- object$mmkin_orig[[1]]$fixed
@@ -200,11 +198,13 @@ print.summary.nlme.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
     print(corr, title = "\nCorrelation:", ...)
   }
 
+  cat("\n") # Random effects
+  print(summary(x$modelStruct), sigma = x$sigma,
+        reEstimates = x$coef$random, verbose = verbose, ...)
+
   cat("\nBacktransformed parameters with asymmetric confidence intervals:\n")
   print(x$confint_back)
 
-  print(summary(x$modelStruct), sigma = x$sigma,
-        reEstimates = x$coef$random, verbose = verbose, ...)
 
   printSFORB <- !is.null(x$SFORB)
   if(printSFORB){
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
new file mode 100644
index 00000000..f7110dd0
--- /dev/null
+++ b/R/summary.saem.mmkin.R
@@ -0,0 +1,250 @@
+#' 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{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 SFO kinetics
+#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+#' dt50_sfo_in_pop <- 50
+#' k_in_pop <- log(2) / dt50_sfo_in_pop
+#' set.seed(1234)
+#' k_in <- rlnorm(5, log(k_in_pop), 0.5)
+#' SFO <- mkinmod(parent = mkinsub("SFO"))
+#'
+#' pred_sfo <- function(k) {
+#'   mkinpredict(SFO,
+#'     c(k_parent = k),
+#'     c(parent = 100),
+#'     sampling_times)
+#' }
+#'
+#' ds_sfo_mean <- lapply(k_in, pred_sfo)
+#' names(ds_sfo_mean) <- paste("ds", 1:5)
+#'
+#' ds_sfo_syn <- lapply(ds_sfo_mean, 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 <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1)
+#' f_saem <- saem(f_mmkin)
+#' summary(f_saem, 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."
+
+  bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
+    object$transform_rates, object$transform_fractions)
+  bpnames <- names(bp)
+
+  #  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 <- object$so@results@name.sigma
+  confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")])
+  colnames(confint_errmod)[1] <- "est."
+
+  # 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
+    }
+  }
+
+  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
+  so_pred <- object$so@results@predictions
+
+  object$data[["observed"]] <- object$data[["value"]]
+  object$data[["value"]] <- NULL
+  object$data[["predicted"]] <- so_pred$ipred
+  object$data[["residual"]] <- so_pred$ires
+  object$data[["standardized"]] <- so_pred$iwres
+  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)
+
+  cat("\nFixed degradation parameter values:\n")
+  if(length(x$fixed$value) == 0) cat("None\n")
+  else print(x$fixed)
+
+  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 = " "))
+
+  cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
+  print(x$confint_trans)
+
+  if (nrow(x$confint_trans) > 1) {
+    corr <- x$corFixed
+    class(corr) <- "correlation"
+    print(corr, title = "\nCorrelation:", ...)
+  }
+
+  cat("\nRandom effects:\n")
+  print(x$confint_ranef)
+
+  cat("\nVariance model:\n")
+  print(x$confint_errmod)
+
+  cat("\nBacktransformed parameters with asymmetric confidence intervals:\n")
+  print(x$confint_back)
+
+  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/_pkgdown.yml b/_pkgdown.yml
index 3668750f..b036ed59 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -45,6 +45,7 @@ reference:
       - plot.nlme.mmkin
       - summary.nlme.mmkin
       - saem.mmkin
+      - summary.saem.mmkin
       - nlme_function
       - get_deg_func
       - saemix_model
diff --git a/man/endpoints.Rd b/man/endpoints.Rd
index 9c354ae9..72487717 100644
--- a/man/endpoints.Rd
+++ b/man/endpoints.Rd
@@ -8,8 +8,11 @@ with mkinfit}
 endpoints(fit)
 }
 \arguments{
-\item{fit}{An object of class \code{\link{mkinfit}} or
-\code{\link{nlme.mmkin}}}
+\item{fit}{An object of class \link{mkinfit}, \link{nlme.mmkin} or
+\link{saem.mmkin}. Or another object that has list components
+mkinmod containing an \link{mkinmod} degradation model, and two numeric vectors,
+bparms.optim and bparms.fixed, that contain parameter values
+for that model.}
 }
 \value{
 A list with a matrix of dissipation times named distimes,
@@ -24,8 +27,13 @@ for one of the parents or metabolites, the Eigenvalues are returned. These
 are equivalent to the rate constants of the DFOP model, but with the
 advantage that the SFORB model can also be used for metabolites.
 }
+\details{
+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
+}
 \note{
-The function is used internally by \code{\link{summary.mkinfit}}.
+The function is used internally by \link{summary.mkinfit},
+\link{summary.nlme.mmkin} and \link{summary.saem.mmkin}.
 }
 \examples{
 
@@ -42,4 +50,3 @@ The function is used internally by \code{\link{summary.mkinfit}}.
 \author{
 Johannes Ranke
 }
-\keyword{manip}
diff --git a/man/saem.Rd b/man/saem.Rd
index 39b66448..96b8b55a 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -43,7 +43,8 @@ used.}
 by the saemix function?}
 }
 \value{
-An \link[saemix:SaemixObject-class]{saemix::SaemixObject}.
+An S3 object of class 'saem.mmkin', containing the fitted
+\link[saemix:SaemixObject-class]{saemix::SaemixObject} as a list component named 'so'.
 
 An \link[saemix:SaemixModel-class]{saemix::SaemixModel} object.
 
@@ -51,8 +52,8 @@ An \link[saemix:SaemixData-class]{saemix::SaemixData} object.
 }
 \description{
 This function uses \code{\link[saemix:saemix]{saemix::saemix()}} as a backend for fitting nonlinear mixed
-effects models created from \link{mmkin} row objects using the stochastic approximation
-to the expectation maximisation algorithm (SAEM).
+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
@@ -95,6 +96,9 @@ f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_typ
 # Using a single core, the following takes about 6 minutes, using 10 cores
 # it is slower instead of faster
 f_saem_des <- saem(f_mmkin_des, cores = 1)
-compare.saemix(list(f_saemix$so, f_saemix_des$so))
+compare.saemix(list(f_saem$so, f_saem_des$so))
 }
 }
+\seealso{
+\link{summary.saem.mmkin}
+}
diff --git a/man/summary.saem.mmkin.Rd b/man/summary.saem.mmkin.Rd
new file mode 100644
index 00000000..0f0d8264
--- /dev/null
+++ b/man/summary.saem.mmkin.Rd
@@ -0,0 +1,89 @@
+% 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}
+\alias{print.summary.saem.mmkin}
+\title{Summary method for class "saem.mmkin"}
+\usage{
+\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{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.saem.mmkin}}
+
+\item{digits}{Number of digits to use for printing}
+}
+\value{
+The summary function returns a list based on the \link[saemix:SaemixObject-class]{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{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 SFO kinetics
+sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+dt50_sfo_in_pop <- 50
+k_in_pop <- log(2) / dt50_sfo_in_pop
+set.seed(1234)
+k_in <- rlnorm(5, log(k_in_pop), 0.5)
+SFO <- mkinmod(parent = mkinsub("SFO"))
+
+pred_sfo <- function(k) {
+  mkinpredict(SFO,
+    c(k_parent = k),
+    c(parent = 100),
+    sampling_times)
+}
+
+ds_sfo_mean <- lapply(k_in, pred_sfo)
+names(ds_sfo_mean) <- paste("ds", 1:5)
+
+ds_sfo_syn <- lapply(ds_sfo_mean, 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 <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1)
+f_saem <- saem(f_mmkin)
+summary(f_saem, data = TRUE)
+}
+
+}
+\author{
+Johannes Ranke for the mkin specific parts
+saemix authors for the parts inherited from saemix.
+}
diff --git a/test.log b/test.log
index 22501b25..32d12209 100644
--- a/test.log
+++ b/test.log
@@ -5,7 +5,7 @@ Testing mkin
 ✔ |   2       | Export dataset for reading into CAKE
 ✔ |  14       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.9 s]
 ✔ |   4       | Calculation of FOCUS chi2 error levels [0.4 s]
-✔ |   7       | Fitting the SFORB model [3.3 s]
+✔ |   7       | Fitting the SFORB model [3.4 s]
 ✔ |   5       | Analytical solutions for coupled models [3.0 s]
 ✔ |   5       | Calculation of Akaike weights
 ✔ |  10       | Confidence intervals and p-values [1.0 s]
@@ -25,7 +25,7 @@ Reason: getRversion() < "4.1.0"  is TRUE
 test_nafta.R:53: skip: Test data from Appendix D are correctly evaluated
 Reason: getRversion() < "4.1.0"  is TRUE
 ────────────────────────────────────────────────────────────────────────────────
-✔ |   9       | Nonlinear mixed-effects models [7.9 s]
+✔ |   9       | Nonlinear mixed-effects models [7.8 s]
 ✔ |   0     1 | Plotting [0.7 s]
 ────────────────────────────────────────────────────────────────────────────────
 test_plot.R:24: skip: Plotting mkinfit and mmkin objects is reproducible
@@ -40,7 +40,7 @@ Reason: getRversion() < "4.1.0"  is TRUE
 ✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.5 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 37.7 s
+Duration: 37.8 s
 
 OK:       146
 Failed:   0
diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf
index dc369ee1..a0e63657 100644
--- a/tests/testthat/FOCUS_2006_D.csf
+++ b/tests/testthat/FOCUS_2006_D.csf
@@ -5,7 +5,7 @@ Description:
 MeasurementUnits: % AR
 TimeUnits: days
 Comments: Created using mkin::CAKE_export
-Date: 2020-11-07
+Date: 2020-11-08
 Optimiser: IRLS
 
 [Data]