Add summary method for nlme.mmkin objects

Improve and update docs

22 files changed, 949 insertions, 93 deletions

diff --git a/NAMESPACE b/NAMESPACE
index 645c145e..9554a0d6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,8 +26,10 @@ S3method(print,mmkin)
 S3method(print,nafta)
 S3method(print,nlme.mmkin)
 S3method(print,summary.mkinfit)
+S3method(print,summary.nlme.mmkin)
 S3method(residuals,mkinfit)
 S3method(summary,mkinfit)
+S3method(summary,nlme.mmkin)
 S3method(update,mkinfit)
 S3method(update,nlme.mmkin)
 export(CAKE_export)
@@ -108,6 +110,7 @@ importFrom(stats,na.fail)
 importFrom(stats,nlminb)
 importFrom(stats,nobs)
 importFrom(stats,optimize)
+importFrom(stats,predict)
 importFrom(stats,pt)
 importFrom(stats,qchisq)
 importFrom(stats,qf)
diff --git a/NEWS.md b/NEWS.md
index ea4b5318..706e0ae8 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,8 +1,6 @@
 # mkin 0.9.50.4 (unreleased)
 
-- 'plot' method for 'nlme.mmkin' objects
-
-- 'print' method for 'mmkin' objects
+- 'plot', 'summary' and 'print' methods for 'nlme.mmkin' objects
 
 - 'saemix_model', 'saemix_data': Helper functions to fit nonlinear mixed-effects models for mmkin row objects using the saemix package
 
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index d3369cf5..6d24a044 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -42,6 +42,9 @@ get_deg_func <- function() {
 #' @importFrom stats na.fail as.formula
 #' @return Upon success, a fitted nlme.mmkin object, which is an nlme object
 #'   with additional elements
+#' @note As the object inherits from [nlme::nlme], there is a wealth of
+#'   methods that will automatically work on 'nlme.mmkin' objects, such as
+#'   [nlme::intervals()], [nlme::anova.lme()] and [nlme::coef.lme()].
 #' @export
 #' @seealso \code{\link{nlme_function}}
 #' @examples
@@ -141,8 +144,8 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
 
   thisCall[["model"]] <- this_model
 
-  mean_dp <- mean_degparms(model)
-  dp_names <- names(mean_dp)
+  mean_dp_start <- mean_degparms(model)
+  dp_names <- names(mean_dp_start)
 
   thisCall[["data"]] <- nlme_data(model)
 
@@ -175,10 +178,21 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
     thisCall[["control"]] <- control
   }
 
-  val <- do.call("nlme.formula", thisCall)
+  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$transform_rates <- model[[1]]$transform_rates
+  val$transform_fractions <- model[[1]]$transform_fractions
+  val$solution_type <- model[[1]]$solution_type
+  val$date.fit <- date()
+  val$nlmeversion <- as.character(utils::packageVersion("nlme"))
+  val$mkinversion <- as.character(utils::packageVersion("mkin"))
+  val$Rversion <- paste(R.version$major, R.version$minor, sep=".")
   class(val) <- c("nlme.mmkin", "nlme", "lme")
   return(val)
 }
@@ -186,10 +200,30 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
 #' @export
 #' @rdname nlme.mmkin
 #' @param x An nlme.mmkin object to print
-#' @param ... Further arguments as in the generic
 print.nlme.mmkin <- function(x, ...) {
-  x$call$data <- "Not shown"
-  NextMethod("print", x)
+  cat( "Kinetic nonlinear mixed-effects model fit by " )
+  cat( if(x$method == "REML") "REML\n" else "maximum likelihood\n")
+  cat("\nStructural model:\n")
+  diffs <- x$mmkin_orig[[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("\nLog-", if(x$method == "REML") "restricted-" else "",
+      "likelihood: ", format(x$logLik), "\n", sep = "")
+  fixF <- x$call$fixed
+  cat("\nFixed effects:\n",
+      deparse(
+  if(inherits(fixF, "formula") || is.call(fixF) || is.name(fixF))
+    x$call$fixed
+  else
+    lapply(fixF, function(el) as.name(deparse(el)))), "\n")
+  print(fixef(x), ...)
+  cat("\n")
+  print(summary(x$modelStruct), sigma = x$sigma, ...)
+  invisible(x)
 }
 
 #' @export
diff --git a/R/plot.nlme.mmkin.R b/R/plot.nlme.mmkin.R
index afb682a7..05a17a22 100644
--- a/R/plot.nlme.mmkin.R
+++ b/R/plot.nlme.mmkin.R
@@ -11,6 +11,8 @@ if(getRversion() >= '2.15.1') utils::globalVariables("ds")
 #' @param rel.height.legend The relative height of the legend shown on top
 #' @param rel.height.bottom The relative height of the bottom plot row
 #' @param ymax Vector of maximum y axis values
+#' @param ncol.legend Number of columns to use in the legend
+#' @param nrow.legend Number of rows to use in the legend
 #' @param \dots Further arguments passed to \code{\link{plot.mkinfit}} and
 #'   \code{\link{mkinresplot}}.
 #' @param resplot Should the residuals plotted against time or against
@@ -28,7 +30,8 @@ if(getRversion() >= '2.15.1') utils::globalVariables("ds")
 #' names(ds) <- paste0("ds ", 6:10)
 #' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
 #'   A1 = mkinsub("SFO"), quiet = TRUE)
-#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, cores = 1)
+#' \dontrun{
+#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
 #' plot(f[, 3:4], standardized = TRUE)
 #'
 #' library(nlme)
@@ -36,6 +39,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables("ds")
 #' # tolerance in order to speed up the fit for this example evaluation
 #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
 #' plot(f_nlme)
+#' }
 #' @export
 plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig),
   obs_vars = names(x$mkinmod$map),
diff --git a/R/sigma_twocomp.R b/R/sigma_twocomp.R
index e8a92ced..e7f4368b 100644
--- a/R/sigma_twocomp.R
+++ b/R/sigma_twocomp.R
@@ -1,16 +1,16 @@
 #' Two-component error model
-#' 
+#'
 #' Function describing the standard deviation of the measurement error in
 #' dependence of the measured value \eqn{y}:
-#' 
+#'
 #' \deqn{\sigma = \sqrt{ \sigma_{low}^2 + y^2 * {rsd}_{high}^2}} sigma =
 #' sqrt(sigma_low^2 + y^2 * rsd_high^2)
-#' 
+#'
 #' This is the error model used for example by Werner et al. (1978). The model
 #' proposed by Rocke and Lorenzato (1995) can be written in this form as well,
 #' but assumes approximate lognormal distribution of errors for high values of
 #' y.
-#' 
+#'
 #' @param y The magnitude of the observed value
 #' @param sigma_low The asymptotic minimum of the standard deviation for low
 #'   observed values
@@ -20,7 +20,7 @@
 #' @references Werner, Mario, Brooks, Samuel H., and Knott, Lancaster B. (1978)
 #'   Additive, Multiplicative, and Mixed Analytical Errors. Clinical Chemistry
 #'   24(11), 1895-1898.
-#' 
+#'
 #'   Rocke, David M. and Lorenzato, Stefan (1995) A two-component model for
 #'   measurement error in analytical chemistry. Technometrics 37(2), 176-184.
 #' @examples
@@ -36,15 +36,8 @@
 #'   data = d_syn, na.action = na.omit,
 #'   start = list(parent_0 = 100, lrc = -3))
 #' if (length(findFunction("varConstProp")) > 0) {
-#'   f_gnls_tc <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
-#'     data = d_syn, na.action = na.omit,
-#'     start = list(parent_0 = 100, lrc = -3),
-#'     weights = varConstProp())
-#'   f_gnls_tc_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
-#'     data = d_syn, na.action = na.omit,
-#'     start = list(parent_0 = 100, lrc = -3),
-#'     control = list(sigma = 1),
-#'     weights = varConstProp())
+#'   f_gnls_tc <- update(f_gnls, weights = varConstProp())
+#'   f_gnls_tc_sf <- update(f_gnls_tc, control = list(sigma = 1))
 #' }
 #' f_mkin <- mkinfit("SFO", d_syn, error_model = "const", quiet = TRUE)
 #' f_mkin_tc <- mkinfit("SFO", d_syn, error_model = "tc", quiet = TRUE)
diff --git a/R/summary.nlme.mmkin.R b/R/summary.nlme.mmkin.R
new file mode 100644
index 00000000..9fdd3f73
--- /dev/null
+++ b/R/summary.nlme.mmkin.R
@@ -0,0 +1,233 @@
+#' Summary method for class "nlme.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 [nlme.mmkin]
+#' @param x an object of class [summary.nlme.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 alpha error level for confidence interval estimation from the t
+#'   distribution
+#' @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 [nlme] object
+#'   obtained in the fit, with at least the following additional components
+#'   \item{nlmeversion, mkinversion, Rversion}{The nlme, 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
+#' @author Johannes Ranke for the mkin specific parts
+#'   José Pinheiro and Douglas Bates for the components inherited from nlme
+#' @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]]
+#' })
+#'
+#' # Evaluate using mmkin and nlme
+#' library(nlme)
+#' f_mmkin <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1)
+#' f_nlme <- nlme(f_mmkin)
+#' summary(f_nlme, data = TRUE)
+#'
+#' @export
+summary.nlme.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, alpha = 0.05, ...) {
+
+  mod_vars <- names(object$mkinmod$diffs)
+
+  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)
+
+  #  variance-covariance estimates for fixed effects (from summary.lme)
+  fixed <- fixef(object)
+  stdFixed <- sqrt(diag(as.matrix(object$varFix)))
+  object$corFixed <- array(
+    t(object$varFix/stdFixed)/stdFixed,
+    dim(object$varFix),
+    list(names(fixed), names(fixed)))
+
+  # 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_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
+  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$verbose <- verbose
+
+  object$fixed <- object$mmkin_orig[[1]]$fixed
+  object$AIC = AIC(object)
+  object$BIC = BIC(object)
+  object$logLik = logLik(object)
+
+  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.nlme.mmkin", "nlme.mmkin", "nlme", "lme")
+  return(object)
+}
+
+#' @rdname summary.nlme.mmkin
+#' @export
+print.summary.nlme.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
+  cat("nlme version used for fitting:     ", x$nlmeversion, "\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", x$numIter, "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")
+  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("\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){
+    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 4631dc97..ac100ac1 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -43,6 +43,7 @@ reference:
     contents:
       - nlme.mmkin
       - plot.nlme.mmkin
+      - summary.nlme.mmkin
       - nlme_function
       - get_deg_func
       - saemix_model
diff --git a/docs/dev/news/index.html b/docs/dev/news/index.html
index c29e413d..4e31f4b7 100644
--- a/docs/dev/news/index.html
+++ b/docs/dev/news/index.html
@@ -146,8 +146,7 @@
 <a href="#mkin-0-9-50-4-unreleased" class="anchor"></a>mkin 0.9.50.4 (unreleased)<small> Unreleased </small>
 </h1>
 <ul>
-<li><p>‘plot’ method for ‘nlme.mmkin’ objects</p></li>
-<li><p>‘print’ method for ‘mmkin’ objects</p></li>
+<li><p>‘plot’, ‘summary’ and ‘print’ methods for ‘nlme.mmkin’ objects</p></li>
 <li><p>‘saemix_model’, ‘saemix_data’: Helper functions to fit nonlinear mixed-effects models for mmkin row objects using the saemix package</p></li>
 </ul>
 </div>
diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml
index 8f493a24..657aa128 100644
--- a/docs/dev/pkgdown.yml
+++ b/docs/dev/pkgdown.yml
@@ -10,7 +10,7 @@ articles:
   web_only/NAFTA_examples: NAFTA_examples.html
   web_only/benchmarks: benchmarks.html
   web_only/compiled_models: compiled_models.html
-last_built: 2020-10-26T13:18Z
+last_built: 2020-10-27T14:34Z
 urls:
   reference: https://pkgdown.jrwb.de/mkin/reference
   article: https://pkgdown.jrwb.de/mkin/articles
diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png
index cfc5bc2b..f001da49 100644
--- a/docs/dev/reference/Rplot001.png
+++ b/docs/dev/reference/Rplot001.png
diff --git a/docs/dev/reference/index.html b/docs/dev/reference/index.html
index 24056025..fa3ec868 100644
--- a/docs/dev/reference/index.html
+++ b/docs/dev/reference/index.html
@@ -308,7 +308,7 @@ of an mmkin object</p></td>
       <tr>
         <th colspan="2">
           <h2 id="section-mixed-models" class="hasAnchor"><a href="#section-mixed-models" class="anchor"></a>Mixed models</h2>
-          <p class="section-desc"><p>Create and work with nonlinear mixed models</p></p>
+          <p class="section-desc"><p>Create and work with nonlinear mixed effects models</p></p>
         </th>
       </tr>
       
@@ -331,6 +331,12 @@ of an mmkin object</p></td>
       </tr><tr>
         
         <td>
+          <p><code><a href="summary.nlme.mmkin.html">summary(<i>&lt;nlme.mmkin&gt;</i>)</a></code> <code><a href="summary.nlme.mmkin.html">print(<i>&lt;summary.nlme.mmkin&gt;</i>)</a></code> </p>
+        </td>
+        <td><p>Summary method for class "nlme.mmkin"</p></td>
+      </tr><tr>
+        
+        <td>
           <p><code><a href="nlme.html">nlme_function()</a></code> <code><a href="nlme.html">mean_degparms()</a></code> <code><a href="nlme.html">nlme_data()</a></code> </p>
         </td>
         <td><p>Helper functions to create nlme models from mmkin row objects</p></td>
diff --git a/docs/dev/reference/nlme.mmkin.html b/docs/dev/reference/nlme.mmkin.html
index 90aec9be..3fa4d97b 100644
--- a/docs/dev/reference/nlme.mmkin.html
+++ b/docs/dev/reference/nlme.mmkin.html
@@ -152,7 +152,7 @@ have been obtained by fitting the same model to a list of datasets.</p>
     </div>
 
     <pre class="usage"><span class='co'># S3 method for mmkin</span>
-<span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span><span class='op'>(</span>
+<span class='fu'>nlme</span><span class='op'>(</span>
   <span class='va'>model</span>,
   data <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/sys.parent.html'>sys.frame</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/sys.parent.html'>sys.parent</a></span><span class='op'>(</span><span class='op'>)</span><span class='op'>)</span>,
   <span class='va'>fixed</span>,
@@ -255,6 +255,11 @@ parameters taken from the mmkin object are used</p></td>
 
     <p>Upon success, a fitted nlme.mmkin object, which is an nlme object
 with additional elements</p>
+    <h2 class="hasAnchor" id="note"><a class="anchor" href="#note"></a>Note</h2>
+
+    <p>As the object inherits from <a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme::nlme</a>, there is a wealth of
+methods that will automatically work on 'nlme.mmkin' objects, such as
+<code><a href='https://rdrr.io/pkg/nlme/man/intervals.html'>nlme::intervals()</a></code>, <code><a href='https://rdrr.io/pkg/nlme/man/anova.lme.html'>nlme::anova.lme()</a></code> and <code><a href='https://rdrr.io/pkg/nlme/man/coef.lme.html'>nlme::coef.lme()</a></code>.</p>
     <h2 class="hasAnchor" id="see-also"><a class="anchor" href="#see-also"></a>See also</h2>
 
     <div class='dont-index'><p><code><a href='nlme_function.html'>nlme_function</a></code></p></div>
@@ -270,11 +275,20 @@ with additional elements</p>
 </div><div class='output co'>#&gt;             df      AIC
 #&gt; f_nlme_sfo   5 625.0539
 #&gt; f_nlme_dfop  9 495.1270</div><div class='input'><span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span><span class='op'>(</span><span class='va'>f_nlme_dfop</span><span class='op'>)</span>
-</div><div class='output co'>#&gt; Nonlinear mixed-effects model fit by maximum likelihood
-#&gt;   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k1,      log_k2, g_ilr) 
-#&gt;   Data: "Not shown" 
-#&gt;   Log-likelihood: -238.5635
-#&gt;   Fixed: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
+</div><div class='output co'>#&gt; Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#&gt; 
+#&gt; Structural model:
+#&gt; d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+#&gt;            time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+#&gt;            * parent
+#&gt; 
+#&gt; Data:
+#&gt; 90 observations of 1 variable(s) grouped in 5 datasets
+#&gt; 
+#&gt; Log-likelihood: -238.5635
+#&gt; 
+#&gt; Fixed effects:
+#&gt;  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
 #&gt;    parent_0      log_k1      log_k2       g_ilr 
 #&gt; 94.17015133 -1.80015306 -4.14738870  0.02290935 
 #&gt; 
@@ -284,9 +298,7 @@ with additional elements</p>
 #&gt;  Structure: Diagonal
 #&gt;         parent_0    log_k1  log_k2     g_ilr Residual
 #&gt; StdDev: 2.488249 0.8447273 1.32965 0.3289311 2.321364
-#&gt; 
-#&gt; Number of Observations: 90
-#&gt; Number of Groups: 5 </div><div class='input'><span class='fu'><a href='https://rdrr.io/r/graphics/plot.default.html'>plot</a></span><span class='op'>(</span><span class='va'>f_nlme_dfop</span><span class='op'>)</span>
+#&gt; </div><div class='input'><span class='fu'><a href='https://rdrr.io/r/graphics/plot.default.html'>plot</a></span><span class='op'>(</span><span class='va'>f_nlme_dfop</span><span class='op'>)</span>
 </div><div class='img'><img src='nlme.mmkin-1.png' alt='' width='700' height='433' /></div><div class='input'><span class='fu'><a href='endpoints.html'>endpoints</a></span><span class='op'>(</span><span class='va'>f_nlme_dfop</span><span class='op'>)</span>
 </div><div class='output co'>#&gt; $distimes
 #&gt;            DT50     DT90 DT50back  DT50_k1  DT50_k2
@@ -295,11 +307,18 @@ with additional elements</p>
 <span class='co'># \dontrun{</span>
   <span class='va'>f_nlme_2</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span><span class='op'>(</span><span class='va'>f</span><span class='op'>[</span><span class='st'>"SFO"</span>, <span class='op'>]</span>, start <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span>parent_0 <span class='op'>=</span> <span class='fl'>100</span>, log_k_parent <span class='op'>=</span> <span class='fl'>0.1</span><span class='op'>)</span><span class='op'>)</span>
   <span class='fu'><a href='https://rdrr.io/r/stats/update.html'>update</a></span><span class='op'>(</span><span class='va'>f_nlme_2</span>, random <span class='op'>=</span> <span class='va'>parent_0</span> <span class='op'>~</span> <span class='fl'>1</span><span class='op'>)</span>
-</div><div class='output co'>#&gt; Nonlinear mixed-effects model fit by maximum likelihood
-#&gt;   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_parent) 
-#&gt;   Data: "Not shown" 
-#&gt;   Log-likelihood: -404.3729
-#&gt;   Fixed: list(parent_0 ~ 1, log_k_parent ~ 1) 
+</div><div class='output co'>#&gt; Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#&gt; 
+#&gt; Structural model:
+#&gt; d_parent/dt = - k_parent * parent
+#&gt; 
+#&gt; Data:
+#&gt;  observations of 0 variable(s) grouped in 0 datasets
+#&gt; 
+#&gt; Log-likelihood: -404.3729
+#&gt; 
+#&gt; Fixed effects:
+#&gt;  list(parent_0 ~ 1, log_k_parent ~ 1) 
 #&gt;     parent_0 log_k_parent 
 #&gt;    75.933480    -3.555983 
 #&gt; 
@@ -307,9 +326,7 @@ with additional elements</p>
 #&gt;  Formula: parent_0 ~ 1 | ds
 #&gt;            parent_0 Residual
 #&gt; StdDev: 0.002416792 21.63027
-#&gt; 
-#&gt; Number of Observations: 90
-#&gt; Number of Groups: 5 </div><div class='input'>  <span class='va'>ds_2</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/lapply.html'>lapply</a></span><span class='op'>(</span><span class='va'>experimental_data_for_UBA_2019</span><span class='op'>[</span><span class='fl'>6</span><span class='op'>:</span><span class='fl'>10</span><span class='op'>]</span>,
+#&gt; </div><div class='input'>  <span class='va'>ds_2</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/lapply.html'>lapply</a></span><span class='op'>(</span><span class='va'>experimental_data_for_UBA_2019</span><span class='op'>[</span><span class='fl'>6</span><span class='op'>:</span><span class='fl'>10</span><span class='op'>]</span>,
    <span class='kw'>function</span><span class='op'>(</span><span class='va'>x</span><span class='op'>)</span> <span class='va'>x</span><span class='op'>$</span><span class='va'>data</span><span class='op'>[</span><span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='st'>"name"</span>, <span class='st'>"time"</span>, <span class='st'>"value"</span><span class='op'>)</span><span class='op'>]</span><span class='op'>)</span>
   <span class='va'>m_sfo_sfo</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span>, <span class='st'>"A1"</span><span class='op'>)</span>,
     A1 <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span><span class='op'>)</span>, use_of_ff <span class='op'>=</span> <span class='st'>"min"</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
@@ -395,11 +412,20 @@ with additional elements</p>
     <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>AIC</a></span><span class='op'>(</span><span class='va'>f_nlme_sfo</span>, <span class='va'>f_nlme_sfo_tc</span>, <span class='va'>f_nlme_dfop</span>, <span class='va'>f_nlme_dfop_tc</span><span class='op'>)</span>
     <span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span><span class='op'>(</span><span class='va'>f_nlme_dfop_tc</span><span class='op'>)</span>
   <span class='op'>}</span>
-</div><div class='output co'>#&gt; Nonlinear mixed-effects model fit by maximum likelihood
-#&gt;   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k1,      log_k2, g_ilr) 
-#&gt;   Data: "Not shown" 
-#&gt;   Log-likelihood: -238.4298
-#&gt;   Fixed: list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
+</div><div class='output co'>#&gt; Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#&gt; 
+#&gt; Structural model:
+#&gt; d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
+#&gt;            time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
+#&gt;            * parent
+#&gt; 
+#&gt; Data:
+#&gt; 90 observations of 1 variable(s) grouped in 5 datasets
+#&gt; 
+#&gt; Log-likelihood: -238.4298
+#&gt; 
+#&gt; Fixed effects:
+#&gt;  list(parent_0 ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_ilr ~ 1) 
 #&gt;    parent_0      log_k1      log_k2       g_ilr 
 #&gt; 94.04774463 -1.82339924 -4.16715509  0.04020161 
 #&gt; 
@@ -415,19 +441,25 @@ with additional elements</p>
 #&gt;  Formula: ~fitted(.) 
 #&gt;  Parameter estimates:
 #&gt;      const       prop 
-#&gt; 2.23222625 0.01262414 
-#&gt; Number of Observations: 90
-#&gt; Number of Groups: 5 </div><div class='input'>
+#&gt; 2.23222625 0.01262414 </div><div class='input'>
   <span class='va'>f_2_obs</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='st'>"SFO-SFO"</span> <span class='op'>=</span> <span class='va'>m_sfo_sfo</span>,
    <span class='st'>"DFOP-SFO"</span> <span class='op'>=</span> <span class='va'>m_dfop_sfo</span><span class='op'>)</span>,
     <span class='va'>ds_2</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span>, error_model <span class='op'>=</span> <span class='st'>"obs"</span><span class='op'>)</span>
   <span class='va'>f_nlme_sfo_sfo_obs</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span><span class='op'>(</span><span class='va'>f_2_obs</span><span class='op'>[</span><span class='st'>"SFO-SFO"</span>, <span class='op'>]</span><span class='op'>)</span>
   <span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span><span class='op'>(</span><span class='va'>f_nlme_sfo_sfo_obs</span><span class='op'>)</span>
-</div><div class='output co'>#&gt; Nonlinear mixed-effects model fit by maximum likelihood
-#&gt;   Model: value ~ (mkin::get_deg_func())(name, time, parent_0, log_k_parent_sink,      log_k_parent_A1, log_k_A1_sink) 
-#&gt;   Data: "Not shown" 
-#&gt;   Log-likelihood: -472.976
-#&gt;   Fixed: list(parent_0 ~ 1, log_k_parent_sink ~ 1, log_k_parent_A1 ~ 1,      log_k_A1_sink ~ 1) 
+</div><div class='output co'>#&gt; Kinetic nonlinear mixed-effects model fit by maximum likelihood
+#&gt; 
+#&gt; Structural model:
+#&gt; d_parent/dt = - k_parent_sink * parent - k_parent_A1 * parent
+#&gt; d_A1/dt = + k_parent_A1 * parent - k_A1_sink * A1
+#&gt; 
+#&gt; Data:
+#&gt; 170 observations of 2 variable(s) grouped in 5 datasets
+#&gt; 
+#&gt; Log-likelihood: -472.976
+#&gt; 
+#&gt; Fixed effects:
+#&gt;  list(parent_0 ~ 1, log_k_parent_sink ~ 1, log_k_parent_A1 ~ 1,      log_k_A1_sink ~ 1) 
 #&gt;          parent_0 log_k_parent_sink   log_k_parent_A1     log_k_A1_sink 
 #&gt;         87.975536         -3.669816         -4.164127         -4.645073 
 #&gt; 
@@ -443,9 +475,7 @@ with additional elements</p>
 #&gt;  Formula: ~1 | name 
 #&gt;  Parameter estimates:
 #&gt;    parent        A1 
-#&gt; 1.0000000 0.2050003 
-#&gt; Number of Observations: 170
-#&gt; Number of Groups: 5 </div><div class='input'>  <span class='co'># The same with DFOP-SFO does not converge, apparently the variances of</span>
+#&gt; 1.0000000 0.2050003 </div><div class='input'>  <span class='co'># The same with DFOP-SFO does not converge, apparently the variances of</span>
   <span class='co'># parent and A1 are too similar in this case, so that the model is</span>
   <span class='co'># overparameterised</span>
   <span class='co'>#f_nlme_dfop_sfo_obs &lt;- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100))</span>
diff --git a/docs/dev/reference/plot.nlme.mmkin.html b/docs/dev/reference/plot.nlme.mmkin.html
index afd9d8d0..267bef05 100644
--- a/docs/dev/reference/plot.nlme.mmkin.html
+++ b/docs/dev/reference/plot.nlme.mmkin.html
@@ -215,6 +215,14 @@ predicted values?</p></td>
 scaling of the y axis and defaults to "auto".</p></td>
     </tr>
     <tr>
+      <th>ncol.legend</th>
+      <td><p>Number of columns to use in the legend</p></td>
+    </tr>
+    <tr>
+      <th>nrow.legend</th>
+      <td><p>Number of rows to use in the legend</p></td>
+    </tr>
+    <tr>
       <th>rel.height.legend</th>
       <td><p>The relative height of the legend shown on top</p></td>
     </tr>
@@ -259,7 +267,8 @@ corresponding model prediction lines for the different datasets.</p></td>
 <span class='fu'><a href='https://rdrr.io/r/base/names.html'>names</a></span><span class='op'>(</span><span class='va'>ds</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/paste.html'>paste0</a></span><span class='op'>(</span><span class='st'>"ds "</span>, <span class='fl'>6</span><span class='op'>:</span><span class='fl'>10</span><span class='op'>)</span>
 <span class='va'>dfop_sfo</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"DFOP"</span>, <span class='st'>"A1"</span><span class='op'>)</span>,
   A1 <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span><span class='op'>)</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
-<span class='va'>f</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='st'>"DFOP-SFO"</span> <span class='op'>=</span> <span class='va'>dfop_sfo</span><span class='op'>)</span>, <span class='va'>ds</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span>, cores <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span>
+<span class='co'># \dontrun{</span>
+<span class='va'>f</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span><span class='st'>"DFOP-SFO"</span> <span class='op'>=</span> <span class='va'>dfop_sfo</span><span class='op'>)</span>, <span class='va'>ds</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
 <span class='fu'><a href='https://rdrr.io/r/graphics/plot.default.html'>plot</a></span><span class='op'>(</span><span class='va'>f</span><span class='op'>[</span>, <span class='fl'>3</span><span class='op'>:</span><span class='fl'>4</span><span class='op'>]</span>, standardized <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
 </div><div class='img'><img src='plot.nlme.mmkin-1.png' alt='' width='700' height='433' /></div><div class='input'>
 <span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'><a href='https://svn.r-project.org/R-packages/trunk/nlme/'>nlme</a></span><span class='op'>)</span>
@@ -267,7 +276,8 @@ corresponding model prediction lines for the different datasets.</p></td>
 <span class='co'># tolerance in order to speed up the fit for this example evaluation</span>
 <span class='va'>f_nlme</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span><span class='op'>(</span><span class='va'>f</span>, control <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>pnlsMaxIter <span class='op'>=</span> <span class='fl'>120</span>, tolerance <span class='op'>=</span> <span class='fl'>1e-3</span><span class='op'>)</span><span class='op'>)</span>
 <span class='fu'><a href='https://rdrr.io/r/graphics/plot.default.html'>plot</a></span><span class='op'>(</span><span class='va'>f_nlme</span><span class='op'>)</span>
-</div><div class='img'><img src='plot.nlme.mmkin-2.png' alt='' width='700' height='433' /></div></pre>
+</div><div class='img'><img src='plot.nlme.mmkin-2.png' alt='' width='700' height='433' /></div><div class='input'><span class='co'># }</span>
+</div></pre>
   </div>
   <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">
     <nav id="toc" data-toggle="toc" class="sticky-top">
diff --git a/docs/dev/reference/sigma_twocomp.html b/docs/dev/reference/sigma_twocomp.html
index ecea8499..8212e480 100644
--- a/docs/dev/reference/sigma_twocomp.html
+++ b/docs/dev/reference/sigma_twocomp.html
@@ -202,15 +202,8 @@ measurement error in analytical chemistry. Technometrics 37(2), 176-184.</p>
   data <span class='op'>=</span> <span class='va'>d_syn</span>, na.action <span class='op'>=</span> <span class='va'>na.omit</span>,
   start <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>parent_0 <span class='op'>=</span> <span class='fl'>100</span>, lrc <span class='op'>=</span> <span class='op'>-</span><span class='fl'>3</span><span class='op'>)</span><span class='op'>)</span>
 <span class='kw'>if</span> <span class='op'>(</span><span class='fu'><a href='https://rdrr.io/r/base/length.html'>length</a></span><span class='op'>(</span><span class='fu'>findFunction</span><span class='op'>(</span><span class='st'>"varConstProp"</span><span class='op'>)</span><span class='op'>)</span> <span class='op'>&gt;</span> <span class='fl'>0</span><span class='op'>)</span> <span class='op'>{</span>
-  <span class='va'>f_gnls_tc</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/gnls.html'>gnls</a></span><span class='op'>(</span><span class='va'>value</span> <span class='op'>~</span> <span class='fu'><a href='https://rdrr.io/r/stats/SSasymp.html'>SSasymp</a></span><span class='op'>(</span><span class='va'>time</span>, <span class='fl'>0</span>, <span class='va'>parent_0</span>, <span class='va'>lrc</span><span class='op'>)</span>,
-    data <span class='op'>=</span> <span class='va'>d_syn</span>, na.action <span class='op'>=</span> <span class='va'>na.omit</span>,
-    start <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>parent_0 <span class='op'>=</span> <span class='fl'>100</span>, lrc <span class='op'>=</span> <span class='op'>-</span><span class='fl'>3</span><span class='op'>)</span>,
-    weights <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/varConstProp.html'>varConstProp</a></span><span class='op'>(</span><span class='op'>)</span><span class='op'>)</span>
-  <span class='va'>f_gnls_tc_sf</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/gnls.html'>gnls</a></span><span class='op'>(</span><span class='va'>value</span> <span class='op'>~</span> <span class='fu'><a href='https://rdrr.io/r/stats/SSasymp.html'>SSasymp</a></span><span class='op'>(</span><span class='va'>time</span>, <span class='fl'>0</span>, <span class='va'>parent_0</span>, <span class='va'>lrc</span><span class='op'>)</span>,
-    data <span class='op'>=</span> <span class='va'>d_syn</span>, na.action <span class='op'>=</span> <span class='va'>na.omit</span>,
-    start <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>parent_0 <span class='op'>=</span> <span class='fl'>100</span>, lrc <span class='op'>=</span> <span class='op'>-</span><span class='fl'>3</span><span class='op'>)</span>,
-    control <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>sigma <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span>,
-    weights <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/varConstProp.html'>varConstProp</a></span><span class='op'>(</span><span class='op'>)</span><span class='op'>)</span>
+  <span class='va'>f_gnls_tc</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/stats/update.html'>update</a></span><span class='op'>(</span><span class='va'>f_gnls</span>, weights <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/varConstProp.html'>varConstProp</a></span><span class='op'>(</span><span class='op'>)</span><span class='op'>)</span>
+  <span class='va'>f_gnls_tc_sf</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/stats/update.html'>update</a></span><span class='op'>(</span><span class='va'>f_gnls_tc</span>, control <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span><span class='op'>(</span>sigma <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span><span class='op'>)</span>
 <span class='op'>}</span>
 <span class='va'>f_mkin</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span><span class='op'>(</span><span class='st'>"SFO"</span>, <span class='va'>d_syn</span>, error_model <span class='op'>=</span> <span class='st'>"const"</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
 <span class='va'>f_mkin_tc</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span><span class='op'>(</span><span class='st'>"SFO"</span>, <span class='va'>d_syn</span>, error_model <span class='op'>=</span> <span class='st'>"tc"</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
diff --git a/docs/dev/reference/summary.nlme.mmkin.html b/docs/dev/reference/summary.nlme.mmkin.html
new file mode 100644
index 00000000..38f53122
--- /dev/null
+++ b/docs/dev/reference/summary.nlme.mmkin.html
@@ -0,0 +1,446 @@
+<!-- Generated by pkgdown: do not edit by hand -->
+<!DOCTYPE html>
+<html lang="en">
+  <head>
+  <meta charset="utf-8">
+<meta http-equiv="X-UA-Compatible" content="IE=edge">
+<meta name="viewport" content="width=device-width, initial-scale=1.0">
+
+<title>Summary method for class "nlme.mmkin" — summary.nlme.mmkin • mkin</title>
+
+
+<!-- jquery -->
+<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.4.1/jquery.min.js" integrity="sha256-CSXorXvZcTkaix6Yvo6HppcZGetbYMGWSFlBw8HfCJo=" crossorigin="anonymous"></script>
+<!-- Bootstrap -->
+
+<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/css/bootstrap.min.css" integrity="sha256-bZLfwXAP04zRMK2BjiO8iu9pf4FbLqX6zitd+tIvLhE=" crossorigin="anonymous" />
+
+<script src="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/js/bootstrap.min.js" integrity="sha256-nuL8/2cJ5NDSSwnKD8VqreErSWHtnEP9E7AySL+1ev4=" crossorigin="anonymous"></script>
+
+<!-- bootstrap-toc -->
+<link rel="stylesheet" href="../bootstrap-toc.css">
+<script src="../bootstrap-toc.js"></script>
+
+<!-- Font Awesome icons -->
+<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/all.min.css" integrity="sha256-mmgLkCYLUQbXn0B1SRqzHar6dCnv9oZFPEC1g1cwlkk=" crossorigin="anonymous" />
+<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/v4-shims.min.css" integrity="sha256-wZjR52fzng1pJHwx4aV2AO3yyTOXrcDW7jBpJtTwVxw=" crossorigin="anonymous" />
+
+<!-- clipboard.js -->
+<script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script>
+
+<!-- headroom.js -->
+<script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/headroom.min.js" integrity="sha256-AsUX4SJE1+yuDu5+mAVzJbuYNPHj/WroHuZ8Ir/CkE0=" crossorigin="anonymous"></script>
+<script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/jQuery.headroom.min.js" integrity="sha256-ZX/yNShbjqsohH1k95liqY9Gd8uOiE1S4vZc+9KQ1K4=" crossorigin="anonymous"></script>
+
+<!-- pkgdown -->
+<link href="../pkgdown.css" rel="stylesheet">
+<script src="../pkgdown.js"></script>
+
+
+
+
+<meta property="og:title" content="Summary method for class "nlme.mmkin" — summary.nlme.mmkin" />
+<meta property="og:description" content="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." />
+
+
+<meta name="robots" content="noindex">
+
+<!-- mathjax -->
+<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js" integrity="sha256-nvJJv9wWKEm88qvoQl9ekL2J+k/RWIsaSScxxlsrv8k=" crossorigin="anonymous"></script>
+<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-84DKXVJXs0/F8OTMzX4UR909+jtl4G7SPypPavF+GfA=" crossorigin="anonymous"></script>
+
+<!--[if lt IE 9]>
+<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
+<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
+<![endif]-->
+
+
+
+  </head>
+
+  <body data-spy="scroll" data-target="#toc">
+    <div class="container template-reference-topic">
+      <header>
+      <div class="navbar navbar-default navbar-fixed-top" role="navigation">
+  <div class="container">
+    <div class="navbar-header">
+      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar" aria-expanded="false">
+        <span class="sr-only">Toggle navigation</span>
+        <span class="icon-bar"></span>
+        <span class="icon-bar"></span>
+        <span class="icon-bar"></span>
+      </button>
+      <span class="navbar-brand">
+        <a class="navbar-link" href="../index.html">mkin</a>
+        <span class="version label label-info" data-toggle="tooltip" data-placement="bottom" title="In-development version">0.9.50.4</span>
+      </span>
+    </div>
+
+    <div id="navbar" class="navbar-collapse collapse">
+      <ul class="nav navbar-nav">
+        <li>
+  <a href="../reference/index.html">Functions and data</a>
+</li>
+<li class="dropdown">
+  <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
+    Articles
+     
+    <span class="caret"></span>
+  </a>
+  <ul class="dropdown-menu" role="menu">
+    <li>
+      <a href="../articles/mkin.html">Introduction to mkin</a>
+    </li>
+    <li>
+      <a href="../articles/FOCUS_D.html">Example evaluation of FOCUS Example Dataset D</a>
+    </li>
+    <li>
+      <a href="../articles/FOCUS_L.html">Example evaluation of FOCUS Laboratory Data L1 to L3</a>
+    </li>
+    <li>
+      <a href="../articles/web_only/FOCUS_Z.html">Example evaluation of FOCUS Example Dataset Z</a>
+    </li>
+    <li>
+      <a href="../articles/web_only/compiled_models.html">Performance benefit by using compiled model definitions in mkin</a>
+    </li>
+    <li>
+      <a href="../articles/twa.html">Calculation of time weighted average concentrations with mkin</a>
+    </li>
+    <li>
+      <a href="../articles/web_only/NAFTA_examples.html">Example evaluation of NAFTA SOP Attachment examples</a>
+    </li>
+    <li>
+      <a href="../articles/web_only/benchmarks.html">Some benchmark timings</a>
+    </li>
+  </ul>
+</li>
+<li>
+  <a href="../news/index.html">News</a>
+</li>
+      </ul>
+      <ul class="nav navbar-nav navbar-right">
+        <li>
+  <a href="https://github.com/jranke/mkin/">
+    <span class="fab fa fab fa-github fa-lg"></span>
+     
+  </a>
+</li>
+      </ul>
+      
+    </div><!--/.nav-collapse -->
+  </div><!--/.container -->
+</div><!--/.navbar -->
+
+      
+
+      </header>
+
+<div class="row">
+  <div class="col-md-9 contents">
+    <div class="page-header">
+    <h1>Summary method for class "nlme.mmkin"</h1>
+    <small class="dont-index">Source: <a href='https://github.com/jranke/mkin/blob/master/R/summary.nlme.mmkin.R'><code>R/summary.nlme.mmkin.R</code></a></small>
+    <div class="hidden name"><code>summary.nlme.mmkin.Rd</code></div>
+    </div>
+
+    <div class="ref-description">
+    <p>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.</p>
+    </div>
+
+    <pre class="usage"><span class='co'># S3 method for nlme.mmkin</span>
+<span class='fu'><a href='https://rdrr.io/r/base/summary.html'>summary</a></span><span class='op'>(</span>
+  <span class='va'>object</span>,
+  data <span class='op'>=</span> <span class='cn'>FALSE</span>,
+  verbose <span class='op'>=</span> <span class='cn'>FALSE</span>,
+  distimes <span class='op'>=</span> <span class='cn'>TRUE</span>,
+  alpha <span class='op'>=</span> <span class='fl'>0.05</span>,
+  <span class='va'>...</span>
+<span class='op'>)</span>
+
+<span class='co'># S3 method for summary.nlme.mmkin</span>
+<span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span><span class='op'>(</span><span class='va'>x</span>, digits <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/Extremes.html'>max</a></span><span class='op'>(</span><span class='fl'>3</span>, <span class='fu'><a href='https://rdrr.io/r/base/options.html'>getOption</a></span><span class='op'>(</span><span class='st'>"digits"</span><span class='op'>)</span> <span class='op'>-</span> <span class='fl'>3</span><span class='op'>)</span>, verbose <span class='op'>=</span> <span class='va'>x</span><span class='op'>$</span><span class='va'>verbose</span>, <span class='va'>...</span><span class='op'>)</span></pre>
+
+    <h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a>Arguments</h2>
+    <table class="ref-arguments">
+    <colgroup><col class="name" /><col class="desc" /></colgroup>
+    <tr>
+      <th>object</th>
+      <td><p>an object of class <a href='nlme.mmkin.html'>nlme.mmkin</a></p></td>
+    </tr>
+    <tr>
+      <th>data</th>
+      <td><p>logical, indicating whether the full data should be included in
+the summary.</p></td>
+    </tr>
+    <tr>
+      <th>verbose</th>
+      <td><p>Should the summary be verbose?</p></td>
+    </tr>
+    <tr>
+      <th>distimes</th>
+      <td><p>logical, indicating whether DT50 and DT90 values should be
+included.</p></td>
+    </tr>
+    <tr>
+      <th>alpha</th>
+      <td><p>error level for confidence interval estimation from the t
+distribution</p></td>
+    </tr>
+    <tr>
+      <th>...</th>
+      <td><p>optional arguments passed to methods like <code>print</code>.</p></td>
+    </tr>
+    <tr>
+      <th>x</th>
+      <td><p>an object of class summary.nlme.mmkin</p></td>
+    </tr>
+    <tr>
+      <th>digits</th>
+      <td><p>Number of digits to use for printing</p></td>
+    </tr>
+    </table>
+
+    <h2 class="hasAnchor" id="value"><a class="anchor" href="#value"></a>Value</h2>
+
+    <p>The summary function returns a list based on the <a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a> object
+obtained in the fit, with at least the following additional components</p>
+<dt>nlmeversion, mkinversion, Rversion</dt><dd><p>The nlme, mkin and R versions used</p></dd>
+<dt>date.fit, date.summary</dt><dd><p>The dates where the fit and the summary were
+produced</p></dd>
+<dt>diffs</dt><dd><p>The differential equations used in the degradation model</p></dd>
+<dt>use_of_ff</dt><dd><p>Was maximum or minimum use made of formation fractions</p></dd>
+<dt>data</dt><dd><p>The data</p></dd>
+<dt>confint_trans</dt><dd><p>Transformed parameters as used in the optimisation, with confidence intervals</p></dd>
+<dt>confint_back</dt><dd><p>Backtransformed parameters, with confidence intervals if available</p></dd>
+<dt>ff</dt><dd><p>The estimated formation fractions derived from the fitted
+model.</p></dd>
+<dt>distimes</dt><dd><p>The DT50 and DT90 values for each observed variable.</p></dd>
+<dt>SFORB</dt><dd><p>If applicable, eigenvalues of SFORB components of the model.</p></dd>
+The print method is called for its side effect, i.e. printing the summary.
+
+    <h2 class="hasAnchor" id="author"><a class="anchor" href="#author"></a>Author</h2>
+
+    <p>Johannes Ranke for the mkin specific parts
+José Pinheiro and Douglas Bates for the components inherited from nlme</p>
+
+    <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
+    <pre class="examples"><div class='input'>
+<span class='co'># Generate five datasets following SFO kinetics</span>
+<span class='va'>sampling_times</span> <span class='op'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span><span class='fl'>0</span>, <span class='fl'>1</span>, <span class='fl'>3</span>, <span class='fl'>7</span>, <span class='fl'>14</span>, <span class='fl'>28</span>, <span class='fl'>60</span>, <span class='fl'>90</span>, <span class='fl'>120</span><span class='op'>)</span>
+<span class='va'>dt50_sfo_in_pop</span> <span class='op'>&lt;-</span> <span class='fl'>50</span>
+<span class='va'>k_in_pop</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/Log.html'>log</a></span><span class='op'>(</span><span class='fl'>2</span><span class='op'>)</span> <span class='op'>/</span> <span class='va'>dt50_sfo_in_pop</span>
+<span class='fu'><a href='https://rdrr.io/r/base/Random.html'>set.seed</a></span><span class='op'>(</span><span class='fl'>1234</span><span class='op'>)</span>
+<span class='va'>k_in</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/stats/Lognormal.html'>rlnorm</a></span><span class='op'>(</span><span class='fl'>5</span>, <span class='fu'><a href='https://rdrr.io/r/base/Log.html'>log</a></span><span class='op'>(</span><span class='va'>k_in_pop</span><span class='op'>)</span>, <span class='fl'>0.5</span><span class='op'>)</span>
+<span class='va'>SFO</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span><span class='op'>(</span><span class='st'>"SFO"</span><span class='op'>)</span><span class='op'>)</span>
+
+<span class='va'>pred_sfo</span> <span class='op'>&lt;-</span> <span class='kw'>function</span><span class='op'>(</span><span class='va'>k</span><span class='op'>)</span> <span class='op'>{</span>
+  <span class='fu'><a href='mkinpredict.html'>mkinpredict</a></span><span class='op'>(</span><span class='va'>SFO</span>,
+    <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span>k_parent <span class='op'>=</span> <span class='va'>k</span><span class='op'>)</span>,
+    <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span><span class='op'>(</span>parent <span class='op'>=</span> <span class='fl'>100</span><span class='op'>)</span>,
+    <span class='va'>sampling_times</span><span class='op'>)</span>
+<span class='op'>}</span>
+
+<span class='va'>ds_sfo_mean</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/lapply.html'>lapply</a></span><span class='op'>(</span><span class='va'>k_in</span>, <span class='va'>pred_sfo</span><span class='op'>)</span>
+<span class='fu'><a href='https://rdrr.io/r/base/names.html'>names</a></span><span class='op'>(</span><span class='va'>ds_sfo_mean</span><span class='op'>)</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/paste.html'>paste</a></span><span class='op'>(</span><span class='st'>"ds"</span>, <span class='fl'>1</span><span class='op'>:</span><span class='fl'>5</span><span class='op'>)</span>
+
+<span class='va'>ds_sfo_syn</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/r/base/lapply.html'>lapply</a></span><span class='op'>(</span><span class='va'>ds_sfo_mean</span>, <span class='kw'>function</span><span class='op'>(</span><span class='va'>ds</span><span class='op'>)</span> <span class='op'>{</span>
+  <span class='fu'><a href='add_err.html'>add_err</a></span><span class='op'>(</span><span class='va'>ds</span>,
+    sdfunc <span class='op'>=</span> <span class='kw'>function</span><span class='op'>(</span><span class='va'>value</span><span class='op'>)</span> <span class='fu'><a href='https://rdrr.io/r/base/MathFun.html'>sqrt</a></span><span class='op'>(</span><span class='fl'>1</span><span class='op'>^</span><span class='fl'>2</span> <span class='op'>+</span> <span class='va'>value</span><span class='op'>^</span><span class='fl'>2</span> <span class='op'>*</span> <span class='fl'>0.07</span><span class='op'>^</span><span class='fl'>2</span><span class='op'>)</span>,
+    n <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span><span class='op'>[[</span><span class='fl'>1</span><span class='op'>]</span><span class='op'>]</span>
+<span class='op'>}</span><span class='op'>)</span>
+
+<span class='co'># Evaluate using mmkin and nlme</span>
+<span class='kw'><a href='https://rdrr.io/r/base/library.html'>library</a></span><span class='op'>(</span><span class='va'><a href='https://svn.r-project.org/R-packages/trunk/nlme/'>nlme</a></span><span class='op'>)</span>
+<span class='va'>f_mmkin</span> <span class='op'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span><span class='op'>(</span><span class='st'>"SFO"</span>, <span class='va'>ds_sfo_syn</span>, quiet <span class='op'>=</span> <span class='cn'>TRUE</span>, error_model <span class='op'>=</span> <span class='st'>"tc"</span>, cores <span class='op'>=</span> <span class='fl'>1</span><span class='op'>)</span>
+</div><div class='output co'>#&gt; <span class='warning'>Warning: Shapiro-Wilk test for standardized residuals: p =  0.00371</span></div><div class='input'><span class='va'>f_nlme</span> <span class='op'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span><span class='op'>(</span><span class='va'>f_mmkin</span><span class='op'>)</span>
+</div><div class='output co'>#&gt; <span class='warning'>Warning: Iteration 3, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)</span></div><div class='input'><span class='fu'><a href='https://rdrr.io/r/base/summary.html'>summary</a></span><span class='op'>(</span><span class='va'>f_nlme</span>, data <span class='op'>=</span> <span class='cn'>TRUE</span><span class='op'>)</span>
+</div><div class='output co'>#&gt; nlme version used for fitting:      3.1.150.1 
+#&gt; mkin version used for pre-fitting:  0.9.50.4 
+#&gt; R version used for fitting:         4.0.3 
+#&gt; Date of fit:     Tue Oct 27 15:25:23 2020 
+#&gt; Date of summary: Tue Oct 27 15:25:23 2020 
+#&gt; 
+#&gt; Equations:
+#&gt; d_parent/dt = - k_parent * parent
+#&gt; 
+#&gt; Data:
+#&gt; 90 observations of 1 variable(s) grouped in 5 datasets
+#&gt; 
+#&gt; Model predictions using solution type analytical 
+#&gt; 
+#&gt; Fitted in 0.581 s using 5 iterations
+#&gt; 
+#&gt; Variance model: Two-component variance function 
+#&gt; 
+#&gt; Mean of starting values for individual parameters:
+#&gt;     parent_0 log_k_parent 
+#&gt;    97.849556    -4.455036 
+#&gt; 
+#&gt; Fixed degradation parameter values:
+#&gt; None
+#&gt; 
+#&gt; Results:
+#&gt; 
+#&gt;       AIC      BIC   logLik
+#&gt;   555.792 570.7908 -271.896
+#&gt; 
+#&gt; Optimised, transformed parameters with symmetric confidence intervals:
+#&gt;                  lower      est.     upper
+#&gt; parent_0     94.701336 97.763446 100.82556
+#&gt; log_k_parent -5.007574 -4.461767  -3.91596
+#&gt; 
+#&gt; Correlation: 
+#&gt;              prnt_0
+#&gt; log_k_parent 0.024 
+#&gt; 
+#&gt; Backtransformed parameters with asymmetric confidence intervals:
+#&gt;                 lower        est.        upper
+#&gt; parent_0 94.701335804 97.76344625 100.82555670
+#&gt; k_parent  0.006687109  0.01154195   0.01992142
+#&gt; Random effects:
+#&gt;  Formula: list(parent_0 ~ 1, log_k_parent ~ 1)
+#&gt;  Level: ds
+#&gt;  Structure: Diagonal
+#&gt;         parent_0 log_k_parent Residual
+#&gt; StdDev: 16.65969     3.516961 5.709013
+#&gt; 
+#&gt; Variance function:
+#&gt;  Structure: Constant plus proportion of variance covariate
+#&gt;  Formula: ~fitted(.) 
+#&gt;  Parameter estimates:
+#&gt;      const       prop 
+#&gt; 1.55075176 0.05680853 
+#&gt; 
+#&gt; Estimated disappearance times:
+#&gt;         DT50  DT90
+#&gt; parent 60.05 199.5
+#&gt; 
+#&gt; Data:
+#&gt;    ds   name time observed predicted residual   std standardized
+#&gt;  ds 1 parent    0    103.6     97.42  6.17540 5.748     1.074413
+#&gt;  ds 1 parent    0     95.9     97.42 -1.52460 5.748    -0.265253
+#&gt;  ds 1 parent    1     95.4     96.72 -1.31719 5.709    -0.230721
+#&gt;  ds 1 parent    1     95.3     96.72 -1.41719 5.709    -0.248237
+#&gt;  ds 1 parent    3     91.6     95.32 -3.71774 5.633    -0.660046
+#&gt;  ds 1 parent    3     94.5     95.32 -0.81774 5.633    -0.145181
+#&gt;  ds 1 parent    7     88.1     92.58 -4.47930 5.483    -0.816920
+#&gt;  ds 1 parent    7     89.6     92.58 -2.97930 5.483    -0.543355
+#&gt;  ds 1 parent   14     90.3     87.97  2.32502 5.233     0.444318
+#&gt;  ds 1 parent   14     96.0     87.97  8.02502 5.233     1.533602
+#&gt;  ds 1 parent   28     80.2     79.44  0.75809 4.772     0.158862
+#&gt;  ds 1 parent   28     77.9     79.44 -1.54191 4.772    -0.323118
+#&gt;  ds 1 parent   60     59.3     62.92 -3.61742 3.896    -0.928458
+#&gt;  ds 1 parent   60     59.6     62.92 -3.31742 3.896    -0.851459
+#&gt;  ds 1 parent   90     59.4     50.56  8.83825 3.264     2.707613
+#&gt;  ds 1 parent   90     51.0     50.56  0.43825 3.264     0.134260
+#&gt;  ds 1 parent  120     38.8     40.63 -1.83247 2.781    -0.658968
+#&gt;  ds 1 parent  120     38.9     40.63 -1.73247 2.781    -0.623007
+#&gt;  ds 2 parent    0    103.2     97.17  6.02995 5.734     1.051655
+#&gt;  ds 2 parent    0     95.1     97.17 -2.07005 5.734    -0.361027
+#&gt;  ds 2 parent    1     88.3     95.59 -7.28901 5.647    -1.290694
+#&gt;  ds 2 parent    1    102.4     95.59  6.81099 5.647     1.206048
+#&gt;  ds 2 parent    3     88.4     92.50 -4.10371 5.479    -0.748984
+#&gt;  ds 2 parent    3     95.2     92.50  2.69629 5.479     0.492110
+#&gt;  ds 2 parent    7     83.5     86.63 -3.12863 5.160    -0.606349
+#&gt;  ds 2 parent    7     96.4     86.63  9.77137 5.160     1.893751
+#&gt;  ds 2 parent   14     77.3     77.23  0.06920 4.653     0.014871
+#&gt;  ds 2 parent   14     76.0     77.23 -1.23080 4.653    -0.264497
+#&gt;  ds 2 parent   28     61.7     61.38  0.31692 3.816     0.083043
+#&gt;  ds 2 parent   28     56.5     61.38 -4.88308 3.816    -1.279513
+#&gt;  ds 2 parent   60     35.1     36.31 -1.21343 2.581    -0.470178
+#&gt;  ds 2 parent   60     32.2     36.31 -4.11343 2.581    -1.593868
+#&gt;  ds 2 parent   90     21.2     22.20 -0.99906 1.999    -0.499832
+#&gt;  ds 2 parent   90     23.3     22.20  1.10094 1.999     0.550800
+#&gt;  ds 2 parent  120     14.1     13.57  0.52931 1.732     0.305638
+#&gt;  ds 2 parent  120     16.9     13.57  3.32931 1.732     1.922443
+#&gt;  ds 3 parent    0     92.4     94.12 -1.71979 5.567    -0.308917
+#&gt;  ds 3 parent    0     94.0     94.12 -0.11979 5.567    -0.021517
+#&gt;  ds 3 parent    1     95.7     91.97  3.72634 5.450     0.683712
+#&gt;  ds 3 parent    1     90.8     91.97 -1.17366 5.450    -0.215343
+#&gt;  ds 3 parent    3     86.7     87.83 -1.12709 5.225    -0.215720
+#&gt;  ds 3 parent    3     85.8     87.83 -2.02709 5.225    -0.387976
+#&gt;  ds 3 parent    7     77.1     80.09 -2.98635 4.807    -0.621300
+#&gt;  ds 3 parent    7     81.5     80.09  1.41365 4.807     0.294104
+#&gt;  ds 3 parent   14     69.1     68.15  0.95467 4.170     0.228922
+#&gt;  ds 3 parent   14     62.4     68.15 -5.74533 4.170    -1.377682
+#&gt;  ds 3 parent   28     49.1     49.34 -0.23911 3.203    -0.074644
+#&gt;  ds 3 parent   28     47.2     49.34 -2.13911 3.203    -0.667787
+#&gt;  ds 3 parent   60     21.9     23.58 -1.68477 2.049    -0.822090
+#&gt;  ds 3 parent   60     23.6     23.58  0.01523 2.049     0.007431
+#&gt;  ds 3 parent   90     12.4     11.81  0.59388 1.690     0.351500
+#&gt;  ds 3 parent   90     13.8     11.81  1.99388 1.690     1.180112
+#&gt;  ds 3 parent  120      4.9      5.91 -1.00993 1.587    -0.636506
+#&gt;  ds 3 parent  120      7.5      5.91  1.59007 1.587     1.002137
+#&gt;  ds 4 parent    0     91.8    101.72 -9.92097 5.983    -1.658171
+#&gt;  ds 4 parent    0    104.6    101.72  2.87903 5.983     0.481194
+#&gt;  ds 4 parent    1    117.5    101.27 16.23017 5.958     2.723944
+#&gt;  ds 4 parent    1     99.3    101.27 -1.96983 5.958    -0.330602
+#&gt;  ds 4 parent    3     94.0    100.37 -6.37355 5.909    -1.078583
+#&gt;  ds 4 parent    3     98.7    100.37 -1.67355 5.909    -0.283212
+#&gt;  ds 4 parent    7    109.2     98.60 10.59529 5.812     1.822915
+#&gt;  ds 4 parent    7     89.2     98.60 -9.40471 5.812    -1.618075
+#&gt;  ds 4 parent   14    103.3     95.58  7.71609 5.647     1.366386
+#&gt;  ds 4 parent   14    103.0     95.58  7.41609 5.647     1.313261
+#&gt;  ds 4 parent   28     90.8     89.82  0.98290 5.333     0.184310
+#&gt;  ds 4 parent   28     88.7     89.82 -1.11710 5.333    -0.209477
+#&gt;  ds 4 parent   60     74.8     77.91 -3.10870 4.690    -0.662879
+#&gt;  ds 4 parent   60     75.3     77.91 -2.60870 4.690    -0.556262
+#&gt;  ds 4 parent   90     71.1     68.18  2.91738 4.172     0.699234
+#&gt;  ds 4 parent   90     78.0     68.18  9.81738 4.172     2.353017
+#&gt;  ds 4 parent  120     59.1     59.67 -0.57073 3.728    -0.153107
+#&gt;  ds 4 parent  120     53.8     59.67 -5.87073 3.728    -1.574902
+#&gt;  ds 5 parent    0     94.9     98.38 -3.48183 5.800    -0.600307
+#&gt;  ds 5 parent    0    101.8     98.38  3.41817 5.800     0.589332
+#&gt;  ds 5 parent    1     96.1     96.75 -0.65141 5.711    -0.114065
+#&gt;  ds 5 parent    1     97.1     96.75  0.34859 5.711     0.061040
+#&gt;  ds 5 parent    3     93.8     93.57  0.22881 5.537     0.041323
+#&gt;  ds 5 parent    3     85.8     93.57 -7.77119 5.537    -1.403444
+#&gt;  ds 5 parent    7     87.6     87.52  0.07909 5.208     0.015186
+#&gt;  ds 5 parent    7     94.0     87.52  6.47909 5.208     1.244026
+#&gt;  ds 5 parent   14     82.5     77.86  4.64101 4.687     0.990182
+#&gt;  ds 5 parent   14     81.7     77.86  3.84101 4.687     0.819498
+#&gt;  ds 5 parent   28     60.0     61.62 -1.61729 3.829    -0.422433
+#&gt;  ds 5 parent   28     61.0     61.62 -0.61729 3.829    -0.161236
+#&gt;  ds 5 parent   60     32.5     36.10 -3.59608 2.571    -1.398750
+#&gt;  ds 5 parent   60     35.5     36.10 -0.59608 2.571    -0.231854
+#&gt;  ds 5 parent   90     21.8     21.86 -0.06415 1.987    -0.032287
+#&gt;  ds 5 parent   90     24.4     21.86  2.53585 1.987     1.276317
+#&gt;  ds 5 parent  120     14.1     13.24  0.85643 1.724     0.496877
+#&gt;  ds 5 parent  120     12.1     13.24 -1.14357 1.724    -0.663473</div><div class='input'>
+</div></pre>
+  </div>
+  <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">
+    <nav id="toc" data-toggle="toc" class="sticky-top">
+      <h2 data-toc-skip>Contents</h2>
+    </nav>
+  </div>
+</div>
+
+
+      <footer>
+      <div class="copyright">
+  <p>Developed by Johannes Ranke.</p>
+</div>
+
+<div class="pkgdown">
+  <p>Site built with <a href="https://pkgdown.r-lib.org/">pkgdown</a> 1.6.1.</p>
+</div>
+
+      </footer>
+   </div>
+
+  
+
+
+  </body>
+</html>
+
+
diff --git a/docs/dev/sitemap.xml b/docs/dev/sitemap.xml
index 8c4d8cbe..67e100de 100644
--- a/docs/dev/sitemap.xml
+++ b/docs/dev/sitemap.xml
@@ -187,6 +187,9 @@
     <loc>https://pkgdown.jrwb.de/mkin/reference/summary.mkinfit.html</loc>
   </url>
   <url>
+    <loc>https://pkgdown.jrwb.de/mkin/reference/summary.nlme.mmkin.html</loc>
+  </url>
+  <url>
     <loc>https://pkgdown.jrwb.de/mkin/reference/synthetic_data_for_UBA_2014.html</loc>
   </url>
   <url>
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
index 041b939a..bf45c6e5 100644
--- a/man/nlme.mmkin.Rd
+++ b/man/nlme.mmkin.Rd
@@ -74,6 +74,11 @@ This functions sets up a nonlinear mixed effects model for an mmkin row
 object. An mmkin row object is essentially a list of mkinfit objects that
 have been obtained by fitting the same model to a list of datasets.
 }
+\note{
+As the object inherits from \link[nlme:nlme]{nlme::nlme}, there is a wealth of
+methods that will automatically work on 'nlme.mmkin' objects, such as
+\code{\link[nlme:intervals]{nlme::intervals()}}, \code{\link[nlme:anova.lme]{nlme::anova.lme()}} and \code{\link[nlme:coef.lme]{nlme::coef.lme()}}.
+}
 \examples{
 ds <- lapply(experimental_data_for_UBA_2019[6:10],
  function(x) subset(x$data[c("name", "time", "value")], name == "parent"))
diff --git a/man/plot.nlme.mmkin.Rd b/man/plot.nlme.mmkin.Rd
index 6944d4b1..5f0f0ef1 100644
--- a/man/plot.nlme.mmkin.Rd
+++ b/man/plot.nlme.mmkin.Rd
@@ -50,6 +50,10 @@ predicted values?}
 \item{maxabs}{Maximum absolute value of the residuals. This is used for the
 scaling of the y axis and defaults to "auto".}
 
+\item{ncol.legend}{Number of columns to use in the legend}
+
+\item{nrow.legend}{Number of rows to use in the legend}
+
 \item{rel.height.legend}{The relative height of the legend shown on top}
 
 \item{rel.height.bottom}{The relative height of the bottom plot row}
@@ -78,7 +82,8 @@ ds <- lapply(experimental_data_for_UBA_2019[6:10],
 names(ds) <- paste0("ds ", 6:10)
 dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
   A1 = mkinsub("SFO"), quiet = TRUE)
-f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, cores = 1)
+\dontrun{
+f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
 plot(f[, 3:4], standardized = TRUE)
 
 library(nlme)
@@ -87,6 +92,7 @@ library(nlme)
 f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
 plot(f_nlme)
 }
+}
 \author{
 Johannes Ranke
 }
diff --git a/man/sigma_twocomp.Rd b/man/sigma_twocomp.Rd
index ed79d493..d205a2f7 100644
--- a/man/sigma_twocomp.Rd
+++ b/man/sigma_twocomp.Rd
@@ -44,15 +44,8 @@ f_gnls <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
   data = d_syn, na.action = na.omit,
   start = list(parent_0 = 100, lrc = -3))
 if (length(findFunction("varConstProp")) > 0) {
-  f_gnls_tc <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
-    data = d_syn, na.action = na.omit,
-    start = list(parent_0 = 100, lrc = -3),
-    weights = varConstProp())
-  f_gnls_tc_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
-    data = d_syn, na.action = na.omit,
-    start = list(parent_0 = 100, lrc = -3),
-    control = list(sigma = 1),
-    weights = varConstProp())
+  f_gnls_tc <- update(f_gnls, weights = varConstProp())
+  f_gnls_tc_sf <- update(f_gnls_tc, control = list(sigma = 1))
 }
 f_mkin <- mkinfit("SFO", d_syn, error_model = "const", quiet = TRUE)
 f_mkin_tc <- mkinfit("SFO", d_syn, error_model = "tc", quiet = TRUE)
diff --git a/man/summary.nlme.mmkin.Rd b/man/summary.nlme.mmkin.Rd
new file mode 100644
index 00000000..ea625dd7
--- /dev/null
+++ b/man/summary.nlme.mmkin.Rd
@@ -0,0 +1,99 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/summary.nlme.mmkin.R
+\name{summary.nlme.mmkin}
+\alias{summary.nlme.mmkin}
+\alias{print.summary.nlme.mmkin}
+\title{Summary method for class "nlme.mmkin"}
+\usage{
+\method{summary}{nlme.mmkin}(
+  object,
+  data = FALSE,
+  verbose = FALSE,
+  distimes = TRUE,
+  alpha = 0.05,
+  ...
+)
+
+\method{print}{summary.nlme.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...)
+}
+\arguments{
+\item{object}{an object of class \link{nlme.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{alpha}{error level for confidence interval estimation from the t
+distribution}
+
+\item{\dots}{optional arguments passed to methods like \code{print}.}
+
+\item{x}{an object of class \link{summary.nlme.mmkin}}
+
+\item{digits}{Number of digits to use for printing}
+}
+\value{
+The summary function returns a list based on the \link{nlme} object
+obtained in the fit, with at least the following additional components
+\item{nlmeversion, mkinversion, Rversion}{The nlme, 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]]
+})
+
+# Evaluate using mmkin and nlme
+library(nlme)
+f_mmkin <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1)
+f_nlme <- nlme(f_mmkin)
+summary(f_nlme, data = TRUE)
+
+}
+\author{
+Johannes Ranke for the mkin specific parts
+José Pinheiro and Douglas Bates for the components inherited from nlme
+}
diff --git a/test.log b/test.log
index 6bc6f08d..10c5645e 100644
--- a/test.log
+++ b/test.log
@@ -5,19 +5,19 @@ 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]
-✔ |   5       | Analytical solutions for coupled models [3.1 s]
+✔ |   7       | Fitting the SFORB model [3.4 s]
+✔ |   5       | Analytical solutions for coupled models [3.2 s]
 ✔ |   5       | Calculation of Akaike weights
 ✔ |  10       | Confidence intervals and p-values [1.0 s]
-✔ |  14       | Error model fitting [3.6 s]
-✔ |   4       | Test fitting the decline of metabolites from their maximum [0.2 s]
+✔ |  14       | Error model fitting [4.2 s]
+✔ |   4       | Test fitting the decline of metabolites from their maximum [0.3 s]
 ✔ |   1       | Fitting the logistic model [0.2 s]
 ✔ |   1       | Test dataset class mkinds used in gmkin
 ✔ |   1       | mkinfit features [0.2 s]
 ✔ |  12       | Special cases of mkinfit calls [0.6 s]
 ✔ |   8       | mkinmod model generation and printing [0.2 s]
 ✔ |   3       | Model predictions with mkinpredict [0.4 s]
-✔ |  14     2 | Evaluations according to 2015 NAFTA guidance [1.0 s]
+✔ |  14     2 | Evaluations according to 2015 NAFTA guidance [1.2 s]
 ────────────────────────────────────────────────────────────────────────────────
 test_nafta.R:25: skip: Test data from Appendix B are correctly evaluated
 Reason: getRversion() < "4.1.0"  is TRUE
@@ -25,8 +25,8 @@ 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.7 s]
-✔ |   0     1 | Plotting [0.6 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
 Reason: getRversion() < "4.1.0"  is TRUE
@@ -35,12 +35,12 @@ Reason: getRversion() < "4.1.0"  is TRUE
 ✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s]
 ✔ |   4       | Summary [0.1 s]
 ✔ |   1       | Summaries of old mkinfit objects
-✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1 s]
-✔ |   9       | Hypothesis tests [6.6 s]
+✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s]
+✔ |   9       | Hypothesis tests [6.7 s]
 ✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.5 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 36.5 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 115d28f4..b8828b57 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-10-08
+Date: 2020-10-27
 Optimiser: IRLS
 
 [Data]