Improve saem method, add summary

Also make the endpoints function work for saem objects.

5 files changed, 341 insertions, 70 deletions

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)
+}