Add summary method for nlme.mmkin objects

Improve and update docs

4 files changed, 285 insertions, 21 deletions

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