mkin coverage - 89.32%

Files
Source

utils::globalVariables(c("predicted", "std"))

#' Fit nonlinear mixed models with SAEM
#'
#' This function uses [saemix::saemix()] as a backend for fitting nonlinear mixed
#' effects models created from [mmkin] row objects using the Stochastic Approximation
#' Expectation Maximisation algorithm (SAEM).
#'
#' An mmkin row object is essentially a list of mkinfit objects that have been
#' obtained by fitting the same model to a list of datasets using [mkinfit].
#'
#' Starting values for the fixed effects (population mean parameters, argument
#' psi0 of [saemix::saemixModel()] are the mean values of the parameters found
#' using [mmkin].
#'
#' @importFrom utils packageVersion
#' @importFrom saemix saemix
#' @param object An [mmkin] row object containing several fits of the same
#' [mkinmod] model to different datasets
#' @param verbose Should we print information about created objects of
#' type [saemix::SaemixModel] and [saemix::SaemixData]?
#' @param transformations Per default, all parameter transformations are done
#' in mkin. If this argument is set to 'saemix', parameter transformations
#' are done in 'saemix' for the supported cases, i.e. (as of version 1.1.2)
#' SFO, FOMC, DFOP and HS without fixing `parent_0`, and SFO or DFOP with
#' one SFO metabolite.
#' @param error_model Possibility to override the error model used in the mmkin object
#' @param degparms_start Parameter values given as a named numeric vector will
#' be used to override the starting values obtained from the 'mmkin' object.
#' @param test_log_parms If TRUE, an attempt is made to use more robust starting
#' values for population parameters fitted as log parameters in mkin (like
#' rate constants) by only considering rate constants that pass the t-test
#' when calculating mean degradation parameters using [mean_degparms].
#' @param conf.level Possibility to adjust the required confidence level
#' for parameter that are tested if requested by 'test_log_parms'.
#' @param solution_type Possibility to specify the solution type in case the
#' automatic choice is not desired
#' @param no_random_effect Character vector of degradation parameters for
#' which there should be no variability over the groups. Only used
#' if the covariance model is not explicitly specified.
#' @param covariance.model Will be passed to [saemix::saemixModel()]. Per
#' default, uncorrelated random effects are specified for all degradation
#' parameters.
#' @param omega.init Will be passed to [saemix::saemixModel()]. If using
#' mkin transformations and the default covariance model with optionally
#' excluded random effects, the variances of the degradation parameters
#' are estimated using [mean_degparms], with testing of untransformed
#' log parameters for significant difference from zero. If not using
#' mkin transformations or a custom covariance model, the default
#' initialisation of [saemix::saemixModel] is used for omega.init.
#' @param covariates A data frame with covariate data for use in
#' 'covariate_models', with dataset names as row names.
#' @param covariate_models A list containing linear model formulas with one explanatory
#' variable, i.e. of the type 'parameter ~ covariate'. Covariates must be available
#' in the 'covariates' data frame.
#' @param error.init Will be passed to [saemix::saemixModel()].
#' @param quiet Should we suppress the messages saemix prints at the beginning
#' and the end of the optimisation process?
#' @param nbiter.saemix Convenience option to increase the number of
#' iterations
#' @param control Passed to [saemix::saemix].
#' @param \dots Further parameters passed to [saemix::saemixModel].
#' @return An S3 object of class 'saem.mmkin', containing the fitted
#' [saemix::SaemixObject] as a list component named 'so'. The
#' object also inherits from 'mixed.mmkin'.
#' @seealso [summary.saem.mmkin] [plot.mixed.mmkin]
#' @examples
#' \dontrun{
#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
#'  function(x) subset(x$data[c("name", "time", "value")]))
#' names(ds) <- paste("Dataset", 6:10)
#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds,
#'   state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
#'
#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop)
#' anova(f_saem_sfo, f_saem_dfop, test = TRUE)
#' illparms(f_saem_dfop)
#' f_saem_dfop_red <- update(f_saem_dfop, no_random_effect = "g_qlogis")
#' anova(f_saem_dfop, f_saem_dfop_red, test = TRUE)
#'
#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop)
#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
#' # functions from saemix
#' library(saemix)
#' compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)
#' plot(f_saem_fomc$so, plot.type = "convergence")
#' plot(f_saem_fomc$so, plot.type = "individual.fit")
#' plot(f_saem_fomc$so, plot.type = "npde")
#' plot(f_saem_fomc$so, plot.type = "vpc")
#'
#' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
#' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
#' anova(f_saem_fomc, f_saem_fomc_tc, test = TRUE)
#'
#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
#'   A1 = mkinsub("SFO"))
#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
#'   A1 = mkinsub("SFO"))
#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
#'   A1 = mkinsub("SFO"))
#' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO,
#' # and compiled ODEs for FOMC that are much slower
#' f_mmkin <- mmkin(list(
#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
#'   ds, quiet = TRUE)
#' # saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds
#' # each on this system, as we use analytical solutions written for saemix.
#' # When using the analytical solutions written for mkin this took around
#' # four minutes
#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
#' f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
#' # We can use print, plot and summary methods to check the results
#' print(f_saem_dfop_sfo)
#' plot(f_saem_dfop_sfo)
#' summary(f_saem_dfop_sfo, data = TRUE)
#'
#' # The following takes about 6 minutes
#' f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
#'   nbiter.saemix = c(200, 80))
#'
#' #anova(
#' #  f_saem_dfop_sfo,
#' #  f_saem_dfop_sfo_deSolve))
#'
#' # If the model supports it, we can also use eigenvalue based solutions, which
#' # take a similar amount of time
#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen",
#' #  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
#' }
#' @export
saem <- function(object, ...) UseMethod("saem")

#' @rdname saem
#' @export
saem.mmkin <- function(object,
  transformations = c("mkin", "saemix"),
  error_model = "auto",
  degparms_start = numeric(),
  test_log_parms = TRUE,
  conf.level = 0.6,
  solution_type = "auto",
  covariance.model = "auto",
  omega.init = "auto",
  covariates = NULL,
  covariate_models = NULL,
  no_random_effect = NULL,
  error.init = c(1, 1),
  nbiter.saemix = c(300, 100),
  control = list(displayProgress = FALSE, print = FALSE,
    nbiter.saemix = nbiter.saemix,
    save = FALSE, save.graphs = FALSE),
  verbose = FALSE, quiet = FALSE, ...)
{
  call <- match.call()
  transformations <- match.arg(transformations)
  m_saemix <- saemix_model(object, verbose = verbose,
    error_model = error_model,
    degparms_start = degparms_start,
    test_log_parms = test_log_parms, conf.level = conf.level,
    solution_type = solution_type,
    transformations = transformations,
    covariance.model = covariance.model,
    omega.init = omega.init,
    covariates = covariates,
    covariate_models = covariate_models,
    error.init = error.init,
    no_random_effect = no_random_effect,
    ...)
  d_saemix <- saemix_data(object, covariates = covariates, verbose = verbose)

  fit_failed <- FALSE
  FIM_failed <- NULL
  fit_time <- system.time({
    utils::capture.output(f_saemix <- try(saemix(m_saemix, d_saemix, control)), split = !quiet)
    if (inherits(f_saemix, "try-error")) fit_failed <- TRUE
  })

  return_data <- nlme_data(object)

  if (!fit_failed) {
    if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects")
    if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) {
      FIM_failed <- c(FIM_failed, "random effects and error model parameters")
    }

    transparms_optim <- f_saemix@results@fixed.effects
    names(transparms_optim) <- f_saemix@results@name.fixed

    if (transformations == "mkin") {
      bparms_optim <- backtransform_odeparms(transparms_optim,
        object[[1]]$mkinmod,
        object[[1]]$transform_rates,
        object[[1]]$transform_fractions)
    } else {
      bparms_optim <- transparms_optim
    }

    saemix_data_ds <- f_saemix@data@data$ds
    mkin_ds_order <- as.character(unique(return_data$ds))
    saemix_ds_order <- unique(saemix_data_ds)

    psi <- saemix::psi(f_saemix)
    rownames(psi) <- saemix_ds_order
    return_data$predicted <- f_saemix@model@model(
      psi = psi[mkin_ds_order, ],
      id = as.numeric(return_data$ds),
      xidep = return_data[c("time", "name")])

    return_data <- transform(return_data,
      residual = value - predicted,
      std = sigma_twocomp(predicted,
        f_saemix@results@respar[1], f_saemix@results@respar[2]))
    return_data <- transform(return_data,
      standardized = residual / std)
  }

  result <- list(
    mkinmod = object[[1]]$mkinmod,
    mmkin = object,
    solution_type = object[[1]]$solution_type,
    transformations = transformations,
    transform_rates = object[[1]]$transform_rates,
    transform_fractions = object[[1]]$transform_fractions,
    covariates = covariates,
    covariate_models = covariate_models,
    sm = m_saemix,
    so = f_saemix,
    call = call,
    time = fit_time,
    FIM_failed = FIM_failed,
    mean_dp_start = attr(m_saemix, "mean_dp_start"),
    bparms.fixed = object[[1]]$bparms.fixed,
    data = return_data,
    err_mod = object[[1]]$err_mod,
    date.fit = date(),
    saemixversion = as.character(utils::packageVersion("saemix")),
    mkinversion = as.character(utils::packageVersion("mkin")),
    Rversion = paste(R.version$major, R.version$minor, sep=".")
  )

  if (!fit_failed) {
    result$mkin_ds_order <- mkin_ds_order
    result$saemix_ds_order <- saemix_ds_order
    result$bparms.optim <- bparms_optim
  }

  class(result) <- c("saem.mmkin", "mixed.mmkin")
  return(result)
}

#' @export
#' @rdname saem
#' @param x An saem.mmkin object to print
#' @param digits Number of digits to use for printing
print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
  cat( "Kinetic nonlinear mixed-effects model fit by SAEM" )
  cat("\nStructural model:\n")
  diffs <- x$mmkin[[1]]$mkinmod$diffs
  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs)
  writeLines(strwrap(nice_diffs, exdent = 11))
  cat("\nData:\n")
  cat(nrow(x$data), "observations of",
    length(unique(x$data$name)), "variable(s) grouped in",
    length(unique(x$data$ds)), "datasets\n")

  if (inherits(x$so, "try-error")) {
    cat("\nFit did not terminate successfully\n")
  } else {
    cat("\nLikelihood computed by importance sampling\n")
    ll <- try(logLik(x$so, type = "is"), silent = TRUE)
    if (inherits(ll, "try-error")) {
      cat("Not available\n")
    } else {
      print(data.frame(
          AIC = AIC(x$so, type = "is"),
          BIC = BIC(x$so, type = "is"),
          logLik = logLik(x$so, type = "is"),
          row.names = " "), digits = digits)
    }

    cat("\nFitted parameters:\n")
    conf.int <- parms(x, ci = TRUE)
    print(conf.int, digits = digits)
  }

  invisible(x)
}

#' @rdname saem
#' @return An [saemix::SaemixModel] object.
#' @export
saemix_model <- function(object, solution_type = "auto",
  transformations = c("mkin", "saemix"), error_model = "auto",
  degparms_start = numeric(),
  covariance.model = "auto", no_random_effect = NULL,
  omega.init = "auto",
  covariates = NULL, covariate_models = NULL,
  error.init = numeric(),
  test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...)
{
  if (nrow(object) > 1) stop("Only row objects allowed")

  mkin_model <- object[[1]]$mkinmod

  if (length(mkin_model$spec) > 1 & solution_type[1] == "analytical") {
    stop("mkin analytical solutions not supported for more thane one observed variable")
  }

  degparms_optim <-  mean_degparms(object, test_log_parms = test_log_parms)
  na_degparms <- names(which(is.na(degparms_optim)))
  if (length(na_degparms) > 0) {
    message("Did not find valid starting values for ", paste(na_degparms, collapse = ", "), "\n",
      "Now trying with test_log_parms = FALSE")
    degparms_optim <-  mean_degparms(object, test_log_parms = FALSE)
  }
  if (transformations == "saemix") {
    degparms_optim <- backtransform_odeparms(degparms_optim,
      object[[1]]$mkinmod,
      object[[1]]$transform_rates,
      object[[1]]$transform_fractions)
  }
  degparms_fixed <- object[[1]]$bparms.fixed

  # Transformations are done in the degradation function by default
  # (transformations = "mkin")
  transform.par = rep(0, length(degparms_optim))

  odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
  odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE)

  odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names)
  odeparms_fixed <- degparms_fixed[odeparms_fixed_names]

  odeini_fixed <- degparms_fixed[odeini_fixed_parm_names]
  names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names)

  model_function <- FALSE

  # Model functions with analytical solutions
  # Fixed parameters, use_of_ff = "min" and turning off sinks currently not supported here
  # In general, we need to consider exactly how the parameters in mkinfit were specified,
  # as the parameters are currently mapped by position in these solutions
  sinks <- sapply(mkin_model$spec, function(x) x$sink)
  if (length(odeparms_fixed) == 0 & mkin_model$use_of_ff == "max" & all(sinks)) {
    # Parent only
    if (length(mkin_model$spec) == 1) {
      parent_type <- mkin_model$spec[[1]]$type
      if (length(odeini_fixed) == 1 && !grepl("_bound$", names(odeini_fixed))) {
        if (transformations == "saemix") {
          stop("saemix transformations are not supported for parent fits with fixed initial parent value")
        }
        if (parent_type == "SFO") {
          stop("saemix needs at least two parameters to work on.")
        }
        if (parent_type == "FOMC") {
          model_function <- function(psi, id, xidep) {
            odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1])
          }
        }
        if (parent_type == "DFOP") {
          model_function <- function(psi, id, xidep) {
            g <- plogis(psi[id, 3])
            t <- xidep[, "time"]
            odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) +
              (1 - g) * exp(- exp(psi[id, 2]) * t))
          }
        }
        if (parent_type == "HS") {
          model_function <- function(psi, id, xidep) {
            tb <- exp(psi[id, 3])
            t <- xidep[, "time"]
            k1 = exp(psi[id, 1])
            odeini_fixed * ifelse(t <= tb,
              exp(- k1 * t),
              exp(- k1 * tb) * exp(- exp(psi[id, 2]) * (t - tb)))
          }
        }
      } else {
        if (length(odeini_fixed) == 2) {
          stop("SFORB with fixed initial parent value is not supported")
        }
        if (parent_type == "SFO") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"])
            }
          } else {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] * exp( - psi[id, 2] * xidep[, "time"])
            }
            transform.par = c(0, 1)
          }
        }
        if (parent_type == "FOMC") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2])
            }
          } else {
            model_function <- function(psi, id, xidep) {
              psi[id, 1] / (xidep[, "time"]/psi[id, 3] + 1)^psi[id, 2]
            }
            transform.par = c(0, 1, 1)
          }
        }
        if (parent_type == "DFOP") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              g <- plogis(psi[id, 4])
              t <- xidep[, "time"]
              psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) +
                (1 - g) * exp(- exp(psi[id, 3]) * t))
            }
          } else {
            model_function <- function(psi, id, xidep) {
              g <- psi[id, 4]
              t <- xidep[, "time"]
              psi[id, 1] * (g * exp(- psi[id, 2] * t) +
                (1 - g) * exp(- psi[id, 3] * t))
            }
            transform.par = c(0, 1, 1, 3)
          }
        }
        if (parent_type == "SFORB") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              k_12 <- exp(psi[id, 3])
              k_21 <- exp(psi[id, 4])
              k_1output <- exp(psi[id, 2])
              t <- xidep[, "time"]

              sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21)
              b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
              b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp

              psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
                ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
            }
          } else {
            model_function <- function(psi, id, xidep) {
              k_12 <- psi[id, 3]
              k_21 <- psi[id, 4]
              k_1output <- psi[id, 2]
              t <- xidep[, "time"]

              sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21)
              b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
              b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp

              psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
                ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
            }
            transform.par = c(0, 1, 1, 1)
          }
        }
        if (parent_type == "HS") {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              tb <- exp(psi[id, 4])
              t <- xidep[, "time"]
              k1 <- exp(psi[id, 2])
              psi[id, 1] * ifelse(t <= tb,
                exp(- k1 * t),
                exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb)))
            }
          } else {
            model_function <- function(psi, id, xidep) {
              tb <- psi[id, 4]
              t <- xidep[, "time"]
              psi[id, 1] * ifelse(t <= tb,
                exp(- psi[id, 2] * t),
                exp(- psi[id, 2] * tb) * exp(- psi[id, 3] * (t - tb)))
            }
            transform.par = c(0, 1, 1, 1)
          }
        }
      }
    }

    # Parent with one metabolite
    # Parameter names used in the model functions are as in
    # https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb
    types <- unname(sapply(mkin_model$spec, function(x) x$type))
    if (length(mkin_model$spec) == 2 &! "SFORB" %in% types ) {
      # Initial value for the metabolite (n20) must be fixed
      if (names(odeini_fixed) == names(mkin_model$spec)[2]) {
        n20 <- odeini_fixed
        parent_name <- names(mkin_model$spec)[1]
        if (identical(types, c("SFO", "SFO"))) {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k1 <- exp(psi[id, 2])
              k2 <- exp(psi[id, 3])
              f12 <- plogis(psi[id, 4])
              ifelse(xidep[, "name"] == parent_name,
                n10 * exp(- k1 * t),
                (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
                  (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
              )
            }
          } else {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k1 <- psi[id, 2]
              k2 <- psi[id, 3]
              f12 <- psi[id, 4]
              ifelse(xidep[, "name"] == parent_name,
                n10 * exp(- k1 * t),
                (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
                  (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
              )
            }
            transform.par = c(0, 1, 1, 3)
          }
        }
        if (identical(types, c("DFOP", "SFO"))) {
          if (transformations == "mkin") {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k2 <- exp(psi[id, 2])
              f12 <- plogis(psi[id, 3])
              l1 <- exp(psi[id, 4])
              l2 <- exp(psi[id, 5])
              g <- plogis(psi[id, 6])
              ifelse(xidep[, "name"] == parent_name,
                n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
                ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
                  (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
                  ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
                      ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
                        f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
                  ((l1 - k2) * l2 - k2 * l1 + k2^2)
              )
            }
          } else {
            model_function <- function(psi, id, xidep) {
              t <- xidep[, "time"]
              n10 <- psi[id, 1]
              k2 <- psi[id, 2]
              f12 <- psi[id, 3]
              l1 <- psi[id, 4]
              l2 <- psi[id, 5]
              g <- psi[id, 6]
              ifelse(xidep[, "name"] == parent_name,
                n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
                ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
                  (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
                  ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
                      ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
                        f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
                  ((l1 - k2) * l2 - k2 * l1 + k2^2)
              )
            }
            transform.par = c(0, 1, 3, 1, 1, 3)
          }
        }
      }
    }
  }

  if (is.function(model_function) & solution_type == "auto") {
    solution_type = "analytical saemix"
  } else {

    if (transformations == "saemix") {
      stop("Using saemix transformations is only supported if an analytical solution is implemented for saemix")
    }

    if (solution_type == "auto")
      solution_type <- object[[1]]$solution_type

    # Define some variables to avoid function calls in model function
    transparms_optim_names <- names(degparms_optim)
    odeini_optim_names <- gsub('_0$', '', odeini_optim_parm_names)
    diff_names <- names(mkin_model$diffs)
    ode_transparms_optim_names <- setdiff(transparms_optim_names, odeini_optim_parm_names)
    transform_rates <- object[[1]]$transform_rates
    transform_fractions <- object[[1]]$transform_fractions

    # Get native symbol info for speed
    use_symbols = FALSE
    if (solution_type == "deSolve" & !is.null(mkin_model$cf)) {
      mkin_model$symbols <- try(deSolve::checkDLL(
        dllname = mkin_model$dll_info[["name"]],
        func = "diffs", initfunc = "initpar",
        jacfunc = NULL, nout = 0, outnames = NULL))
      if (!inherits(mkin_model$symbols, "try-error")) {
        use_symbols = TRUE
      }
    }

    # Define the model function
    model_function <- function(psi, id, xidep) {

      uid <- unique(id)

      res_list <- lapply(uid, function(i) {

        transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict
        names(transparms_optim) <- transparms_optim_names

        odeini_optim <- transparms_optim[odeini_optim_parm_names]
        names(odeini_optim) <- odeini_optim_names

        odeini <- c(odeini_optim, odeini_fixed)[diff_names]

        odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model,
          transform_rates = transform_rates,
          transform_fractions = transform_fractions)
        odeparms <- c(odeparms_optim, odeparms_fixed)

        xidep_i <- xidep[which(id == i), ]

        if (solution_type[1] == "analytical") {
          out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms)
        } else {

          i_time <- xidep_i$time
          i_name <- xidep_i$name

          out_wide <- mkinpredict(mkin_model,
            odeparms = odeparms, odeini = odeini,
            solution_type = solution_type,
            outtimes = sort(unique(i_time)),
            na_stop = FALSE
          )

          out_index <- cbind(as.character(i_time), as.character(i_name))
          out_values <- out_wide[out_index]
        }
        return(out_values)
      })
      res <- unlist(res_list)
      return(res)
    }
  }

  if (identical(error_model, "auto")) {
    error_model = object[[1]]$err_mod
  }
  error.model <- switch(error_model,
    const = "constant",
    tc = "combined",
    obs = "constant")

  if (error_model == "obs") {
    warning("The error model 'obs' (variance by variable) can currently not be transferred to an saemix model")
  }

  degparms_psi0 <- degparms_optim
  degparms_psi0[names(degparms_start)] <- degparms_start
  psi0_matrix <- matrix(degparms_psi0, nrow = 1,
    dimnames = list("(Intercept)", names(degparms_psi0)))

  if (covariance.model[1] == "auto") {
    covariance_diagonal <- rep(1, length(degparms_optim))
    if (!is.null(no_random_effect)) {
      degparms_no_random <- which(names(degparms_psi0) %in% no_random_effect)
      covariance_diagonal[degparms_no_random] <- 0
    }
    covariance.model = diag(covariance_diagonal)
  }

  if (omega.init[1] == "auto") {
    if (transformations == "mkin") {
      degparms_eta_ini <- as.numeric( # remove names
        mean_degparms(object,
          random = TRUE, test_log_parms = TRUE)$eta)

      omega.init <- 2 * diag(degparms_eta_ini^2)
    } else {
      omega.init <- matrix(nrow = 0, ncol = 0)
    }
  }

  if (is.null(covariate_models)) {
    covariate.model <- matrix(nrow = 0, ncol = 0) # default in saemixModel()
  } else {
    degparms_dependent <- sapply(covariate_models, function(x) as.character(x[[2]]))
    covariates_in_models = unique(unlist(lapply(
      covariate_models, function(x)
        colnames(attr(terms(x), "factors"))
      )))
    covariates_not_available <- setdiff(covariates_in_models, names(covariates))
    if (length(covariates_not_available) > 0) {
      stop("Covariate(s) ", paste(covariates_not_available, collapse = ", "),
        " used in the covariate models not available in the covariate data")
    }
    psi0_matrix <- rbind(psi0_matrix,
      matrix(0, nrow = length(covariates), ncol = ncol(psi0_matrix),
        dimnames = list(names(covariates), colnames(psi0_matrix))))
    covariate.model <- matrix(0, nrow = length(covariates),
      ncol = ncol(psi0_matrix),
      dimnames = list(
        covariates = names(covariates),
        degparms = colnames(psi0_matrix)))
    if (transformations == "saemix") {
      stop("Covariate models with saemix transformations currently not supported")
    }
    parms_trans <- as.data.frame(t(sapply(object, parms, transformed = TRUE)))
    for (covariate_model in covariate_models) {
      covariate_name <- as.character(covariate_model[[2]])
      model_data <- cbind(parms_trans, covariates)
      ini_model <- lm(covariate_model, data = model_data)
      ini_coef <- coef(ini_model)
      psi0_matrix[names(ini_coef), covariate_name] <- ini_coef
      covariate.model[names(ini_coef)[-1], covariate_name] <- as.numeric(as.logical(ini_coef[-1]))
    }
  }

  res <- saemix::saemixModel(model_function,
    psi0 = psi0_matrix,
    "Mixed model generated from mmkin object",
    transform.par = transform.par,
    error.model = error.model,
    verbose = verbose,
    covariance.model = covariance.model,
    covariate.model = covariate.model,
    omega.init = omega.init,
    error.init = error.init,
    ...
  )
  attr(res, "mean_dp_start") <- degparms_optim
  return(res)
}

#' @rdname saem
#' @importFrom rlang !!!
#' @return An [saemix::SaemixData] object.
#' @export
saemix_data <- function(object, covariates = NULL, verbose = FALSE, ...) {
  if (nrow(object) > 1) stop("Only row objects allowed")
  ds_names <- colnames(object)

  ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
  names(ds_list) <- ds_names
  ds_saemix_all <- vctrs::vec_rbind(!!!ds_list, .names_to = "ds")
  ds_saemix <- data.frame(ds = ds_saemix_all$ds,
    name = as.character(ds_saemix_all$variable),
    time = ds_saemix_all$time,
    value = ds_saemix_all$observed,
    stringsAsFactors = FALSE)
  if (!is.null(covariates)) {
    name.covariates <- names(covariates)
    covariates$ds <- rownames(covariates)
    ds_saemix <- merge(ds_saemix, covariates, sort = FALSE)
  } else {
    name.covariates <- character(0)
  }

  res <- saemix::saemixData(ds_saemix,
    name.group = "ds",
    name.predictors = c("time", "name"),
    name.response = "value",
    name.covariates = name.covariates,
    verbose = verbose,
    ...)
  return(res)
}

#' logLik method for saem.mmkin objects
#'
#' @param object The fitted [saem.mmkin] object
#' @param \dots Passed to [saemix::logLik.SaemixObject]
#' @param method Passed to [saemix::logLik.SaemixObject]
#' @export
logLik.saem.mmkin <- function(object, ..., method = c("is", "lin", "gq")) {
  method <- match.arg(method)
  return(logLik(object$so, method = method))
}

#' @export
update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
  call <- object$call
  # For some reason we get saem.mmkin in the call when using mhmkin
  # so we need to fix this so we do not have to export saem.mmkin in
  # addition to the S3 method
  call[[1]] <- saem

  # We also need to provide the mmkin object in the call, so it
  # will also be found when called by testthat or pkgdown
  call[[2]] <- object$mmkin

  update_arguments <- match.call(expand.dots = FALSE)$...

  if (length(update_arguments) > 0) {
    update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
  }

  for (a in names(update_arguments)[update_arguments_in_call]) {
    call[[a]] <- update_arguments[[a]]
  }

  update_arguments_not_in_call <- !update_arguments_in_call
  if(any(update_arguments_not_in_call)) {
    call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
    call <- as.call(call)
  }
  if(evaluate) eval(call, parent.frame())
  else call
}

#' @export
#' @rdname parms
#' @param ci Should a matrix with estimates and confidence interval boundaries
#' be returned? If FALSE (default), a vector of estimates is returned if no
#' covariates are given, otherwise a matrix of estimates is returned, with
#' each column corresponding to a row of the data frame holding the covariates
#' @param covariates A data frame holding covariate values for which to
#' return parameter values. Only has an effect if 'ci' is FALSE.
parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, ...) {
  cov.mod <- object$sm@covariance.model
  n_cov_mod_parms <- sum(cov.mod[upper.tri(cov.mod, diag = TRUE)])
  n_parms <- length(object$sm@name.modpar) +
    n_cov_mod_parms +
    length(object$sm@name.sigma)

  if (inherits(object$so, "try-error")) {
    conf.int <- matrix(rep(NA, 3 * n_parms), ncol = 3)
    colnames(conf.int) <- c("estimate", "lower", "upper")
  } else {
    conf.int <- object$so@results@conf.int[c("estimate", "lower", "upper")]
    rownames(conf.int) <- object$so@results@conf.int[["name"]]
    conf.int.var <- grepl("^Var\\.", rownames(conf.int))
    conf.int <- conf.int[!conf.int.var, ]
    conf.int.cov <- grepl("^Cov\\.", rownames(conf.int))
    conf.int <- conf.int[!conf.int.cov, ]
  }
  estimate <- conf.int[, "estimate"]

  names(estimate) <- rownames(conf.int)

  if (ci) {
    return(conf.int)
  } else {
    if (is.null(covariates)) {
      return(estimate)
    } else {
      est_for_cov <- matrix(NA,
        nrow = length(object$sm@name.modpar), ncol = nrow(covariates),
        dimnames = (list(object$sm@name.modpar, rownames(covariates))))
      covmods <- object$covariate_models
      names(covmods) <- sapply(covmods, function(x) as.character(x[[2]]))
      for (deg_parm_name in rownames(est_for_cov)) {
        if (deg_parm_name %in% names(covmods)) {
          covariate <- covmods[[deg_parm_name]][[3]]
          beta_degparm_name <- paste0("beta_", covariate,
            "(", deg_parm_name, ")")
          est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name] +
            covariates[[covariate]] * estimate[beta_degparm_name]
        } else {
          est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name]
        }
      }
      return(est_for_cov)
    }
  }
}

#' Method to get status information for fit array objects
#'
#' @param object The object to investigate
#' @param x The object to be printed
#' @param \dots For potential future extensions
#' @return  An object with the same dimensions as the fit array
#' suitable printing method.
#' @export
status <- function(object, ...)
{
  UseMethod("status", object)
}

#' @rdname status
#' @export
#' @examples
#' \dontrun{
#' fits <- mmkin(
#'   c("SFO", "FOMC"),
#'   list("FOCUS A" = FOCUS_2006_A,
#'        "FOCUS B" = FOCUS_2006_C),
#'   quiet = TRUE)
#' status(fits)
#' }
status.mmkin <- function(object, ...) {
  all_summary_warnings <- character()
  sww <- 0 # Counter for Shapiro-Wilks warnings

  result <- lapply(object,
    function(fit) {
      if (inherits(fit, "try-error")) return("E")
      sw <- fit$summary_warnings
      swn <- names(sw)
      if (length(sw) > 0) {
        if (any(grepl("S", swn))) {
          sww <<- sww + 1
          swn <- gsub("S", paste0("S", sww), swn)
        }
        warnstring <- paste(swn, collapse = ", ")
        names(sw) <- swn
        all_summary_warnings <<- c(all_summary_warnings, sw)
        return(warnstring)
      } else {
        return("OK")
      }
    })
  result <- unlist(result)
  dim(result) <- dim(object)
  dimnames(result) <- dimnames(object)

  u_swn <- unique(names(all_summary_warnings))
  attr(result, "unique_warnings") <- all_summary_warnings[u_swn]
  class(result) <- "status.mmkin"
  return(result)
}

#' @rdname status
#' @export
print.status.mmkin <- function(x, ...) {
  u_w <- attr(x, "unique_warnings")
  attr(x, "unique_warnings") <- NULL
  class(x) <- NULL
  print(x, quote = FALSE)
  cat("\n")
  for (i in seq_along(u_w)) {
    cat(names(u_w)[i], ": ", u_w[i], "\n", sep = "")
  }
  if (any(x == "OK")) cat("OK: No warnings\n")
  if (any(x == "E")) cat("E: Error\n")
}

#' @rdname status
#' @export
status.mhmkin <- function(object, ...) {
  if (inherits(object[[1]], "saem.mmkin")) {
    test_func <- function(fit) {
      if (inherits(fit, "try-error")) {
          return("E")
      } else {
        if (inherits(fit$so, "try-error")) {
          return("E")
        } else {
          if (!is.null(fit$FIM_failed)) {
            return_values <- c("fixed effects" = "Fth",
              "random effects and error model parameters" = "FO")
            return(paste(return_values[fit$FIM_failed], collapse = ", "))
          } else {
            return("OK")
          }
        }
      }
    }
  } else {
    stop("Only mhmkin objects containing saem.mmkin objects currently supported")
  }
  result <- lapply(object, test_func)
  result <- unlist(result)
  dim(result) <- dim(object)
  dimnames(result) <- dimnames(object)

  class(result) <- "status.mhmkin"
  return(result)
}

#' @rdname status
#' @export
print.status.mhmkin <- function(x, ...) {
  class(x) <- NULL
  print(x, quote = FALSE)
  cat("\n")
  if (any(x == "OK")) cat("OK: Fit terminated successfully\n")
  if (any(x == "Fth")) cat("Fth: Could not invert FIM for fixed effects\n")
  if (any(x == "FO")) cat("FO: Could not invert FIM for random effects and error model parameters\n")
  if (any(x == "Fth, FO")) cat("Fth, FO: Could not invert FIM for fixed effects, nor for random effects and error model parameters\n")
  if (any(x == "E")) cat("E: Error\n")
}


#' Function to calculate endpoints for further use from kinetic models fitted
#' with mkinfit
#'
#' This function calculates DT50 and DT90 values as well as formation fractions
#' from kinetic models fitted with mkinfit. If the SFORB model was specified
#' for one of the parents or metabolites, the Eigenvalues are returned. These
#' are equivalent to the rate constants of the DFOP model, but with the
#' advantage that the SFORB model can also be used for metabolites.
#'
#' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
#' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
#'
#' @param fit An object of class [mkinfit], [nlme.mmkin] or [saem.mmkin], or
#' another object that has list components mkinmod containing an [mkinmod]
#' degradation model, and two numeric vectors, bparms.optim and bparms.fixed,
#' that contain parameter values for that model.
#' @param covariates Numeric vector with covariate values for all variables in
#' any covariate models in the object. If given, it overrides 'covariate_quantile'.
#' @param covariate_quantile This argument only has an effect if the fitted
#' object has covariate models. If so, the default is to show endpoints
#' for the median of the covariate values (50th percentile).
#' @importFrom stats optimize
#' @return A list with a matrix of dissipation times named distimes, and, if
#' applicable, a vector of formation fractions named ff and, if the SFORB model
#' was in use, a vector of eigenvalues of these SFORB models, equivalent to
#' DFOP rate constants
#' @note The function is used internally by [summary.mkinfit],
#' [summary.nlme.mmkin] and [summary.saem.mmkin].
#' @author Johannes Ranke
#' @examples
#'
#'   fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
#'   endpoints(fit)
#'   \dontrun{
#'     fit_2 <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE)
#'     endpoints(fit_2)
#'     fit_3 <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE)
#'     endpoints(fit_3)
#'   }
#'
#' @export
endpoints <- function(fit, covariates = NULL, covariate_quantile = 0.5) {
  mkinmod <- fit$mkinmod
  obs_vars <- names(mkinmod$spec)

  if (!is.null(fit$covariate_models)) {
    if (is.null(covariates)) {
      covariates = as.data.frame(
        apply(fit$covariates, 2, quantile,
          covariate_quantile, simplify = FALSE))
    } else {
      covariate_m <- matrix(covariates, byrow = TRUE)
      colnames(covariate_m) <- names(covariates)
      rownames(covariate_m) <- "User"
      covariates <- as.data.frame(covariate_m)
    }
    degparms_trans <- parms(fit, covariates = covariates)[, 1]
    if (inherits(fit, "saem.mmkin") & (fit$transformations == "saemix")) {
      degparms <- degparms_trans
    } else {
      degparms <- backtransform_odeparms(degparms_trans,
        fit$mkinmod,
        transform_rates = fit$transform_rates,
        transform_fractions = fit$transform_fractions)
    }
  } else {
    degparms <- c(fit$bparms.optim, fit$bparms.fixed)
  }

  # Set up object to return
  ep <- list()
  ep$covariates <- covariates
  ep$ff <- vector()
  ep$SFORB <- vector()
  ep$distimes <- data.frame(
    DT50 = rep(NA, length(obs_vars)),
    DT90 = rep(NA, length(obs_vars)),
    row.names = obs_vars)

  for (obs_var in obs_vars) {
    type = names(mkinmod$map[[obs_var]])[1]

    # Get formation fractions if directly fitted, and calculate remaining fraction to sink
    f_names = grep(paste("^f", obs_var, sep = "_"), names(degparms), value=TRUE)
    if (length(f_names) > 0) {
      f_values = degparms[f_names]
      f_to_sink = 1 - sum(f_values)
      names(f_to_sink) = ifelse(type == "SFORB",
        paste(obs_var, "free", "sink", sep = "_"),
        paste(obs_var, "sink", sep = "_"))
      for (f_name in f_names) {
        ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
      }
      ep$ff = append(ep$ff, f_to_sink)
    }

    # Get the rest
    if (type == "SFO") {
      k_names = grep(paste("^k", obs_var, sep="_"), names(degparms), value=TRUE)
      k_tot = sum(degparms[k_names])
      DT50 = log(2)/k_tot
      DT90 = log(10)/k_tot
      if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) {
        for (k_name in k_names)
        {
          ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
        }
      }
    }
    if (type == "FOMC") {
      alpha = degparms["alpha"]
      beta = degparms["beta"]
      DT50 = beta * (2^(1/alpha) - 1)
      DT90 = beta * (10^(1/alpha) - 1)
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
      ep$distimes[obs_var, c("DT50back")] = DT50_back
    }
    if (type == "IORE") {
      k_names = grep(paste("^k__iore", obs_var, sep="_"), names(degparms), value=TRUE)
      k_tot = sum(degparms[k_names])
      # From the NAFTA kinetics guidance, p. 5
      n = degparms[paste("N", obs_var, sep = "_")]
      k = k_tot
      # Use the initial concentration of the parent compound
      source_name = mkinmod$map[[1]][[1]]
      c0 = degparms[paste(source_name, "0", sep = "_")]
      alpha = 1 / (n - 1)
      beta = (c0^(1 - n))/(k * (n - 1))
      DT50 = beta * (2^(1/alpha) - 1)
      DT90 = beta * (10^(1/alpha) - 1)
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
      ep$distimes[obs_var, c("DT50back")] = DT50_back
      if (mkinmod$use_of_ff == "min") {
        for (k_name in k_names)
        {
          ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
        }
      }
    }
    if (type == "DFOP") {
      k1 = degparms["k1"]
      k2 = degparms["k2"]
      g = degparms["g"]
      f <- function(log_t, x) {
        t <- exp(log_t)
        fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
        (fraction - (1 - x/100))^2
      }
      DT50_k1 = log(2)/k1
      DT50_k2 = log(2)/k2
      DT90_k1 = log(10)/k1
      DT90_k2 = log(10)/k2

      DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum),
                  silent = TRUE)
      DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum),
                  silent = TRUE)
      if (inherits(DT50, "try-error")) DT50 = NA
      if (inherits(DT90, "try-error")) DT90 = NA
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011

      ep$distimes[obs_var, c("DT50back")] = DT50_back
      ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
      ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
    }
    if (type == "HS") {
      k1 = degparms["k1"]
      k2 = degparms["k2"]
      tb = degparms["tb"]
      DTx <- function(x) {
        DTx.a <- (log(100/(100 - x)))/k1
        DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
        if (DTx.a < tb) DTx <- DTx.a
        else DTx <- DTx.b
        return(DTx)
      }
      DT50 <- DTx(50)
      DT90 <- DTx(90)
      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
      DT50_k1 = log(2)/k1
      DT50_k2 = log(2)/k2
      ep$distimes[obs_var, c("DT50back")] = DT50_back
      ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
      ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
    }
    if (type == "SFORB") {
      # FOCUS kinetics (2006), p. 60 f
      k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(degparms), value=TRUE)
      k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
      k_1output = sum(degparms[k_out_names])
      k_12 = degparms[paste("k", obs_var, "free", "bound", sep="_")]
      k_21 = degparms[paste("k", obs_var, "bound", "free", sep="_")]

      sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21)
      b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
      b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
      g = (k_12 + k_21 - b1)/(b2 - b1)

      DT50_b1 = log(2)/b1
      DT50_b2 = log(2)/b2
      DT90_b1 = log(10)/b1
      DT90_b2 = log(10)/b2

      SFORB_fraction = function(t) {
        g * exp(-b1 * t) + (1 - g) * exp(-b2 * t)
      }

      f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2
      log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum,
                      silent = TRUE)
      f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2
      log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum,
                      silent = TRUE)

      DT50 = if (inherits(log_DT50, "try-error")) NA
             else exp(log_DT50)
      DT90 = if (inherits(log_DT90, "try-error")) NA
             else exp(log_DT90)

      DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011

      for (k_out_name in k_out_names)
      {
        ep$ff[[sub("k_", "", k_out_name)]] = degparms[[k_out_name]] / k_1output
      }

      # Return the eigenvalues for comparison with DFOP rate constants
      ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
      ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
      # Return g for comparison with DFOP
      ep$SFORB[[paste(obs_var, "g", sep="_")]] = g

      ep$distimes[obs_var, c("DT50back")] = DT50_back
      ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1
      ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2
    }
    if (type == "logistic") {
      # FOCUS kinetics (2014) p. 67
      kmax = degparms["kmax"]
      k0 = degparms["k0"]
      r = degparms["r"]
      DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))
      DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))

      DT50_k0 = log(2)/k0
      DT50_kmax = log(2)/kmax
      ep$distimes[obs_var, c("DT50_k0")] = DT50_k0
      ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax
    }
    ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)
  }
  if (length(ep$ff) == 0) ep$ff <- NULL
  if (length(ep$SFORB) == 0) ep$SFORB <- NULL
  return(ep)
}

#' Fit nonlinear mixed-effects models built from one or more kinetic
#' degradation models and one or more error models
#'
#' The name of the methods expresses that (**m**ultiple) **h**ierarchichal
#' (also known as multilevel) **m**ulticompartment **kin**etic models are
#' fitted. Our kinetic models are nonlinear, so we can use various nonlinear
#' mixed-effects model fitting functions.
#'
#' @param objects A list of [mmkin] objects containing fits of the same
#' degradation models to the same data, but using different error models.
#' Alternatively, a single [mmkin] object containing fits of several
#' degradation models to the same data
#' @param backend The backend to be used for fitting. Currently, only saemix is
#' supported
#' @param no_random_effect Default is NULL and will be passed to [saem]. If a
#' character vector is supplied, it will be passed to all calls to [saem],
#' which will exclude random effects for all matching parameters. Alternatively,
#' a list of character vectors or an object of class [illparms.mhmkin] can be
#' specified. They have to have the same dimensions that the return object of
#' the current call will have, i.e. the number of rows must match the number
#' of degradation models in the mmkin object(s), and the number of columns must
#' match the number of error models used in the mmkin object(s).
#' @param algorithm The algorithm to be used for fitting (currently not used)
#' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
#' model fitting function.
#' @param cores The number of cores to be used for multicore processing. This
#' is only used when the \code{cluster} argument is \code{NULL}. On Windows
#' machines, cores > 1 is not supported, you need to use the \code{cluster}
#' argument to use multiple logical processors. Per default, all cores detected
#' by [parallel::detectCores()] are used, except on Windows where the default
#' is 1.
#' @param cluster A cluster as returned by [makeCluster] to be used for
#' parallel execution.
#' @importFrom parallel mclapply parLapply detectCores
#' @return A two-dimensional [array] of fit objects and/or try-errors that can
#' be indexed using the degradation model names for the first index (row index)
#' and the error model names for the second index (column index), with class
#' attribute 'mhmkin'.
#' @author Johannes Ranke
#' @seealso \code{\link{[.mhmkin}} for subsetting [mhmkin] objects
#' @export
mhmkin <- function(objects, ...) {
  UseMethod("mhmkin")
}

#' @export
#' @rdname mhmkin
mhmkin.mmkin <- function(objects, ...) {
  mhmkin(list(objects), ...)
}

#' @export
#' @rdname mhmkin
#' @examples
#' \dontrun{
#' # We start with separate evaluations of all the first six datasets with two
#' # degradation models and two error models
#' f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
#' f_sep_tc <- update(f_sep_const, error_model = "tc")
#' # The mhmkin function sets up hierarchical degradation models aka
#' # nonlinear mixed-effects models for all four combinations, specifying
#' # uncorrelated random effects for all degradation parameters
#' f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
#' status(f_saem_1)
#' # The 'illparms' function shows that in all hierarchical fits, at least
#' # one random effect is ill-defined (the confidence interval for the
#' # random effect expressed as standard deviation includes zero)
#' illparms(f_saem_1)
#' # Therefore we repeat the fits, excluding the ill-defined random effects
#' f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
#' status(f_saem_2)
#' illparms(f_saem_2)
#' # Model comparisons show that FOMC with two-component error is preferable,
#' # and confirms our reduction of the default parameter model
#' anova(f_saem_1)
#' anova(f_saem_2)
#' # The convergence plot for the selected model looks fine
#' saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
#' # The plot of predictions versus data shows that we have a pretty data-rich
#' # situation with homogeneous distribution of residuals, because we used the
#' # same degradation model, error model and parameter distribution model that
#' # was used in the data generation.
#' plot(f_saem_2[["FOMC", "tc"]])
#' # We can specify the same parameter model reductions manually
#' no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
#' dim(no_ranef) <- c(2, 2)
#' f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
#' anova(f_saem_2m)
#' }
mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
  no_random_effect = NULL,
  ...,
  cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
{
  call <- match.call()
  dot_args <- list(...)
  backend_function <- switch(backend,
    saemix = "saem"
  )

  deg_models <- lapply(objects[[1]][, 1], function(x) x$mkinmod)
  names(deg_models) <- dimnames(objects[[1]])$model
  n.deg <- length(deg_models)

  ds <- lapply(objects[[1]][1, ], function(x) x$data)

  for (other in objects[-1]) {
    # Check if the degradation models in all objects are the same
    for (deg_model_name in names(deg_models)) {
      if (!all.equal(other[[deg_model_name, 1]]$mkinmod$spec,
            deg_models[[deg_model_name]]$spec))
      {
        stop("The mmkin objects have to be based on the same degradation models")
      }
    }
    # Check if they have been fitted to the same dataset
    other_object_ds <- lapply(other[1, ], function(x) x$data)
    for (i in 1:length(ds)) {
      if (!all.equal(ds[[i]][c("time", "variable", "observed")],
          other_object_ds[[i]][c("time", "variable", "observed")]))
      {
        stop("The mmkin objects have to be fitted to the same datasets")
      }
    }
  }

  n.o <- length(objects)

  error_models = sapply(objects, function(x) x[[1]]$err_mod)
  n.e <- length(error_models)

  n.fits <- n.deg * n.e
  fit_indices <- matrix(1:n.fits, ncol = n.e)
  dimnames(fit_indices) <- list(degradation = names(deg_models),
                                error = error_models)

  if (is.null(no_random_effect) || is.null(dim(no_random_effect))) {
    no_ranef <- rep(list(no_random_effect), n.fits)
    dim(no_ranef) <- dim(fit_indices)
  } else {
    if (!identical(dim(no_random_effect), dim(fit_indices))) {
      stop("Dimensions of argument 'no_random_effect' are not suitable")
    }
    if (is(no_random_effect, "illparms.mhmkin")) {
      no_ranef_dim <- dim(no_random_effect)
      no_ranef <- lapply(no_random_effect, function(x) {
        no_ranef_split <- strsplit(x, ", ")
        ret <- sapply(no_ranef_split, function(y) {
          gsub("sd\\((.*)\\)", "\\1", y)
        })
        return(ret)
      })
      dim(no_ranef) <- no_ranef_dim
    } else {
      no_ranef <- no_random_effect
    }
  }

  fit_function <- function(fit_index) {
    w <- which(fit_indices == fit_index, arr.ind = TRUE)
    deg_index <- w[1]
    error_index <- w[2]
    mmkin_row <- objects[[error_index]][deg_index, ]
    res <- try(do.call(backend_function,
        args = c(
          list(mmkin_row),
          dot_args,
          list(no_random_effect = no_ranef[[deg_index, error_index]]))))
    return(res)
  }


  fit_time <- system.time({
    if (is.null(cluster)) {
      results <- parallel::mclapply(as.list(1:n.fits), fit_function,
        mc.cores = cores, mc.preschedule = FALSE)
    } else {
      results <- parallel::parLapply(cluster, as.list(1:n.fits), fit_function)
    }
  })

  attributes(results) <- attributes(fit_indices)
  attr(results, "call") <- call
  attr(results, "time") <- fit_time
  class(results) <- switch(backend,
    saemix = c("mhmkin.saem.mmkin", "mhmkin")
  )
  return(results)
}

#' Subsetting method for mhmkin objects
#'
#' @param x An [mhmkin] object.
#' @param i Row index selecting the fits for specific models
#' @param j Column index selecting the fits to specific datasets
#' @param drop If FALSE, the method always returns an mhmkin object, otherwise
#'   either a list of fit objects or a single fit object.
#' @return An object inheriting from \code{\link{mhmkin}}.
#' @rdname mhmkin
#' @export
`[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
  original_class <- class(x)
  class(x) <- NULL
  x_sub <- x[i, j, drop = drop]

  if (!drop) {
    class(x_sub) <- original_class
  }
  return(x_sub)
}

#' Print method for mhmkin objects
#'
#' @rdname mhmkin
#' @export
print.mhmkin <- function(x, ...) {
  cat("<mhmkin> object\n")
  cat("Status of individual fits:\n\n")
  print(status(x))
}

#' Check if fit within an mhmkin object failed
#' @param x The object to be checked
check_failed <- function(x) {
  if (inherits(x, "try-error")) {
    return(TRUE)
  } else {
    if (inherits(x$so, "try-error")) {
      return(TRUE)
    } else {
      return(FALSE)
    }
  }
}

#' @export
AIC.mhmkin <- function(object, ..., k = 2) {
  res <- sapply(object, function(x) {
    if (check_failed(x)) return(NA)
    else return(AIC(x$so, k = k))
  })
  dim(res) <- dim(object)
  dimnames(res) <- dimnames(object)
  return(res)
}

#' @export
BIC.mhmkin <- function(object, ...) {
  res <- sapply(object, function(x) {
    if (check_failed(x)) return(NA)
    else return(BIC(x$so))
  })
  dim(res) <- dim(object)
  dimnames(res) <- dimnames(object)
  return(res)
}

#' @export
update.mhmkin <- function(object, ..., evaluate = TRUE) {
  call <- attr(object, "call")
  # For some reason we get mhkin.list in call[[1]] when using mhmkin from the
  # loaded package so we need to fix this so we do not have to export
  # mhmkin.list in addition to the S3 method mhmkin
  call[[1]] <- mhmkin

  update_arguments <- match.call(expand.dots = FALSE)$...

  if (length(update_arguments) > 0) {
    update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
  }

  for (a in names(update_arguments)[update_arguments_in_call]) {
    call[[a]] <- update_arguments[[a]]
  }

  update_arguments_not_in_call <- !update_arguments_in_call
  if(any(update_arguments_not_in_call)) {
    call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
    call <- as.call(call)
  }
  if(evaluate) eval(call, parent.frame())
  else call
}

#' @export
anova.mhmkin <- function(object, ...,
  method = c("is", "lin", "gq"), test = FALSE, model.names = "auto") {
  if (identical(model.names, "auto")) {
    model.names <- outer(rownames(object), colnames(object), paste)
  }
  failed_index <- which(sapply(object, check_failed), arr.ind = TRUE)
  if (length(failed_index > 0)) {
    rlang::inject(anova(!!!(object[-failed_index]), method = method, test = test,
      model.names = model.names[-failed_index]))
  } else {
    rlang::inject(anova(!!!(object), method = method, test = test,
      model.names = model.names))
  }
}


#' Evaluate parent kinetics using the NAFTA guidance
#'
#' The function fits the SFO, IORE and DFOP models using \code{\link{mmkin}}
#' and returns an object of class \code{nafta} that has methods for printing
#' and plotting.
#'
#' @param ds A dataframe that must contain one variable called "time" with the
#'   time values specified by the \code{time} argument, one column called
#'   "name" with the grouping of the observed values, and finally one column of
#'   observed values called "value".
#' @param title Optional title of the dataset
#' @param quiet Should the evaluation text be shown?
#' @param \dots Further arguments passed to \code{\link{mmkin}} (not for the
#'  printing method).
#' @importFrom stats qf
#' @return An list of class \code{nafta}. The list element named "mmkin" is the
#'   \code{\link{mmkin}} object containing the fits of the three models. The
#'   list element named "title" contains the title of the dataset used. The
#'   list element "data" contains the dataset used in the fits.
#' @author Johannes Ranke
#' @source NAFTA (2011) Guidance for evaluating and calculating degradation
#'   kinetics in environmental media. NAFTA Technical Working Group on
#'   Pesticides
#'   \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/guidance-evaluating-and-calculating-degradation}
#'   accessed 2019-02-22
#'
#'   US EPA (2015) Standard Operating Procedure for Using the NAFTA Guidance to
#'   Calculate Representative Half-life Values and Characterizing Pesticide
#'   Degradation
#'   \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/standard-operating-procedure-using-nafta-guidance}
#' @examples
#'
#'   nafta_evaluation <- nafta(NAFTA_SOP_Appendix_D, cores = 1)
#'   print(nafta_evaluation)
#'   plot(nafta_evaluation)
#'
#' @export
nafta <- function(ds, title = NA, quiet = FALSE, ...) {
  if (length(levels(ds$name)) > 1) {
    stop("The NAFTA procedure is only defined for decline data for a single compound")
  }
  n <- nrow(subset(ds, !is.na(value)))
  models <- c("SFO", "IORE", "DFOP")

  result <- list(title = title, data = ds)
  result$mmkin <- mmkin(models, list(ds), quiet = TRUE, ...)

  distimes <- lapply(result$mmkin, function(x) as.numeric(endpoints(x)$distimes["parent", ]))

  result$distimes <- matrix(NA, nrow = 3, ncol = 3,
    dimnames = list(models, c("DT50", "DT90", "DT50_rep")))
  result$distimes["SFO", ] <- distimes[[1]][c(1, 2, 1)]
  result$distimes["IORE", ] <- distimes[[2]][c(1, 2, 3)]
  result$distimes["DFOP", ] <- distimes[[3]][c(1, 2, 5)]

  # Get parameters with statistics
  result$parameters <- lapply(result$mmkin, function(x) {
    summary(x)$bpar[, c(1, 4:6)]
  })
  names(result$parameters) <- models

  # Compare the sum of squared residuals (SSR) to the upper bound of the
  # confidence region of the SSR for the IORE model
  result$S <- sapply(result$mmkin, function(x) sum(x$data$residual^2))
  names(result$S) <- c("SFO", "IORE", "DFOP")
  # Equation (3) on p. 3
  p <- 3
  result$S["IORE"]
  result$S_c <- result$S[["IORE"]] * (1 + p/(n - p) * qf(0.5, p, n - p))

  result$t_rep <- .evaluate_nafta_results(result$S, result$S_c,
    result$distimes, quiet = quiet)

  class(result) <- "nafta"
  return(result)
}

#' Plot the results of the three models used in the NAFTA scheme.
#'
#' The plots are ordered with increasing complexity of the model in this
#' function (SFO, then IORE, then DFOP).
#'
#' Calls \code{\link{plot.mmkin}}.
#'
#' @param x An object of class \code{\link{nafta}}.
#' @param legend Should a legend be added?
#' @param main Possibility to override the main title of the plot.
#' @param \dots Further arguments passed to \code{\link{plot.mmkin}}.
#' @return The function is called for its side effect.
#' @author Johannes Ranke
#' @export
plot.nafta <- function(x, legend = FALSE, main = "auto", ...) {
  if (main == "auto") {
    if (is.na(x$title)) main = ""
    else main = x$title
  }
  plot(x$mmkin, ..., legend = legend, main = main)
}

#' Print nafta objects
#'
#' Print nafta objects. The results for the three models are printed in the
#' order of increasing model complexity, i.e. SFO, then IORE, and finally DFOP.
#'
#' @param x An \code{\link{nafta}} object.
#' @param digits Number of digits to be used for printing parameters and
#'   dissipation times.
#' @rdname nafta
#' @export
print.nafta <- function(x, quiet = TRUE, digits = 3, ...) {
  cat("Sums of squares:\n")
  print(x$S)
  cat("\nCritical sum of squares for checking the SFO model:\n")
  print(x$S_c)
  cat("\nParameters:\n")
  print(x$parameters, digits = digits)
  t_rep <- .evaluate_nafta_results(x$S, x$S_c, x$distimes, quiet = quiet)
  cat("\nDTx values:\n")
  print(signif(x$distimes, digits = digits))
  cat("\nRepresentative half-life:\n")
  print(round(t_rep, 2))
}

.evaluate_nafta_results <- function(S, S_c, distimes, quiet = FALSE) {
  t_SFO <- distimes["IORE", "DT50"]
  t_IORE <- distimes["IORE", "DT50_rep"]
  t_DFOP2 <- distimes["DFOP", "DT50_rep"]

  if (S["SFO"] < S_c) {
    if (!quiet) {
      message("S_SFO is lower than the critical value S_c, use the SFO model")
    }
    t_rep <- t_SFO
  } else {
    if (!quiet) {
      message("The SFO model is rejected as S_SFO is equal or higher than the critical value S_c")
    }
    if (t_IORE < t_DFOP2) {
      if (!quiet) {
        message("The half-life obtained from the IORE model may be used")
      }
      t_rep <- t_IORE
    } else {
      if (!quiet) {
        message("The representative half-life of the IORE model is longer than the one corresponding")
        message("to the terminal degradation rate found with the DFOP model.")
        message("The representative half-life obtained from the DFOP model may be used")
      }
      t_rep <- t_DFOP2
    }
  }
  return(t_rep)
}

#' 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)
#'
#' set.seed(12345)
#' ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) {
#'   add_err(ds,
#'     sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
#'     n = 1)[[1]]
#' })
#'
#' \dontrun{
#' # Evaluate using mmkin and 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)

  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[[1]]$error_model_algorithm
  err_mod = object$mmkin[[1]]$err_mod

  object$diffs <- object$mkinmod$diffs
  object$print_data <- data
  object$data[["observed"]] <- object$data[["value"]]
  object$data[["value"]] <- NULL
  object$data[["predicted"]] <- predict(object)
  object$data[["residual"]] <- residuals(object, type = "response")
  if (is.null(object$modelStruct$varStruct)) {
    object$data[["std"]] <- object$sigma
  } else {
    object$data[["std"]] <- 1/attr(object$modelStruct$varStruct, "weights")
  }
  object$data[["standardized"]] <- residuals(object, type = "pearson")
  object$verbose <- verbose

  object$fixed <- object$mmkin[[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, digits = digits)

  cat("\nFixed degradation parameter values:\n")
  if(length(x$fixed$value) == 0) cat("None\n")
  else print(x$fixed, digits = digits)

  cat("\nResults:\n\n")
  print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
      row.names = " "), digits = digits, ...)

  cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
  print(x$confint_trans, digits = digits, ...)

  if (nrow(x$confint_trans) > 1) {
    corr <- x$corFixed
    class(corr) <- "correlation"
    print(corr, title = "\nCorrelation:", rdig = digits, ...)
  }

  cat("\n") # Random effects
  print(summary(x$modelStruct), sigma = x$sigma,
        reEstimates = x$coef$random, digits = digits, verbose = verbose, ...)

  cat("\nBacktransformed parameters with asymmetric confidence intervals:\n")
  print(x$confint_back, digits = digits, ...)


  printSFORB <- !is.null(x$SFORB)
  if(printSFORB){
    cat("\nEstimated Eigenvalues of SFORB model(s):\n")
    print(x$SFORB, digits = digits,...)
  }

  printff <- !is.null(x$ff)
  if(printff){
    cat("\nResulting formation fractions:\n")
    print(data.frame(ff = x$ff), digits = digits, ...)
  }

  printdistimes <- !is.null(x$distimes)
  if(printdistimes){
    cat("\nEstimated disappearance times:\n")
    print(x$distimes, digits = digits, ...)
  }

  if (x$print_data){
    cat("\nData:\n")
    print(format(x$data, digits = digits, ...), row.names = FALSE)
  }

  invisible(x)
}

utils::globalVariables(c("variable", "residual"))

#' Function to plot squared residuals and the error model for an mkin object
#' 
#' This function plots the squared residuals for the specified subset of the
#' observed variables from an mkinfit object. In addition, one or more dashed
#' line(s) show the fitted error model.  A combined plot of the fitted model
#' and this error model plot can be obtained with \code{\link{plot.mkinfit}}
#' using the argument \code{show_errplot = TRUE}.
#' 
#' @param object A fit represented in an \code{\link{mkinfit}} object.
#' @param obs_vars A character vector of names of the observed variables for
#'   which residuals should be plotted. Defaults to all observed variables in
#'   the model
#' @param xlim plot range in x direction.
#' @param xlab Label for the x axis.
#' @param ylab Label for the y axis.
#' @param maxy Maximum value of the residuals. This is used for the scaling of
#'   the y axis and defaults to "auto".
#' @param legend Should a legend be plotted?
#' @param lpos Where should the legend be placed? Default is "topright". Will
#'   be passed on to \code{\link{legend}}.
#' @param col_obs Colors for the observed variables.
#' @param pch_obs Symbols to be used for the observed variables.
#' @param frame Should a frame be drawn around the plots?
#' @param \dots further arguments passed to \code{\link{plot}}.
#' @return Nothing is returned by this function, as it is called for its side
#'   effect, namely to produce a plot.
#' @author Johannes Ranke
#' @seealso \code{\link{mkinplot}}, for a way to plot the data and the fitted
#'   lines of the mkinfit object.
#' @keywords hplot
#' @examples
#' 
#' \dontrun{
#' model <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
#' fit <- mkinfit(model, FOCUS_2006_D, error_model = "tc", quiet = TRUE)
#' mkinerrplot(fit)
#' }
#' 
#' @export
mkinerrplot <- function (object,
  obs_vars = names(object$mkinmod$map),
  xlim = c(0, 1.1 * max(object$data$predicted)),
  xlab = "Predicted", ylab = "Squared residual",
  maxy = "auto", legend= TRUE, lpos = "topright",
  col_obs = "auto", pch_obs = "auto",
  frame = TRUE,
  ...)
{
  obs_vars_all <- as.character(unique(object$data$variable))

  if (length(obs_vars) > 0){
      obs_vars <- intersect(obs_vars_all, obs_vars)
  } else obs_vars <- obs_vars_all

  residuals <- subset(object$data, variable %in% obs_vars, residual)

  if (maxy == "auto") maxy = max(residuals^2, na.rm = TRUE)

  # Set colors and symbols
  if (col_obs[1] == "auto") {
    col_obs <- 1:length(obs_vars)
  }

  if (pch_obs[1] == "auto") {
    pch_obs <- 1:length(obs_vars)
  }
  names(col_obs) <- names(pch_obs) <- obs_vars

  plot(0, type = "n",
       xlab = xlab, ylab = ylab,
       xlim = xlim,
       ylim = c(0, 1.2 * maxy), frame = frame, ...)

  for(obs_var in obs_vars){
    residuals_plot <- subset(object$data, variable == obs_var, c("predicted", "residual"))
    points(residuals_plot[["predicted"]],
           residuals_plot[["residual"]]^2,
           pch = pch_obs[obs_var], col = col_obs[obs_var])
  }

  if (object$err_mod == "const") {
    abline(h = object$errparms^2, lty = 2, col = 1)
  }
  if (object$err_mod == "obs") {
    for (obs_var in obs_vars) {
      sigma_name = paste0("sigma_", obs_var)
      abline(h = object$errparms[sigma_name]^2, lty = 2,
             col = col_obs[obs_var])
    }
  }
  if (object$err_mod == "tc") {
    sigma_plot <- function(predicted) {
      sigma_twocomp(predicted,
                    sigma_low = object$errparms[1],
                    rsd_high = object$errparms[2])^2
    }
    plot(sigma_plot, from = 0, to = max(object$data$predicted),
         add = TRUE, lty = 2, col = 1)
  }

  if (legend == TRUE) {
    legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
      col = col_obs[obs_vars], pch = pch_obs[obs_vars])
  }
}

#' Plot parameter variability of multistart objects
#'
#' Produces a boxplot with all parameters from the multiple runs, scaled
#' either by the parameters of the run with the highest likelihood,
#' or by their medians as proposed in the paper by Duchesne et al. (2021).
#'
#' Starting values of degradation model parameters and error model parameters
#' are shown as green circles. The results obtained in the original run
#' are shown as red circles.
#'
#' @param object The [multistart] object
#' @param llmin The minimum likelihood of objects to be shown
#' @param llquant Fractional value for selecting only the fits with higher
#' likelihoods. Overrides 'llmin'.
#' @param scale By default, scale parameters using the best
#' available fit.
#' If 'median', parameters are scaled using the median parameters from all fits.
#' @param main Title of the plot
#' @param lpos Positioning of the legend.
#' @param \dots Passed to [boxplot]
#' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical
#' identifiability in the frame of nonlinear mixed effects models: the example
#' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478.
#' doi: 10.1186/s12859-021-04373-4.
#' @seealso [multistart]
#' @importFrom stats median quantile
#' @export
parplot <- function(object, ...) {
  UseMethod("parplot")
}

#' @rdname parplot
#' @export
parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
  scale = c("best", "median"),
  lpos = "bottomleft", main = "", ...)
{
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar, no.readonly = TRUE))

  orig <- attr(object, "orig")
  orig_parms <- parms(orig)
  start_degparms <- orig$mean_dp_start
  all_parms <- parms(object, exclude_failed = FALSE)

  if (inherits(object, "multistart.saem.mmkin")) {
    llfunc <- function(object) {
      if (inherits(object$so, "try-error")) return(NA)
      else return(logLik(object$so))
    }
  } else {
    stop("parplot is only implemented for multistart.saem.mmkin objects")
  }
  ll <- sapply(object, llfunc)
  if (!is.na(llquant[1])) {
    if (llmin != -Inf) warning("Overriding 'llmin' because 'llquant' was specified")
    llmin <- quantile(ll, 1 - llquant)
  }
  selected <- which(ll > llmin)
  selected_parms <- all_parms[selected, ]

  par(las = 1)
  if (orig$transformations == "mkin") {
    degparm_names_transformed <- names(start_degparms)
    degparm_index <- which(names(orig_parms) %in% degparm_names_transformed)
    orig_parms[degparm_names_transformed] <- backtransform_odeparms(
      orig_parms[degparm_names_transformed],
      orig$mmkin[[1]]$mkinmod,
      transform_rates = orig$mmkin[[1]]$transform_rates,
      transform_fractions = orig$mmkin[[1]]$transform_fractions)
    start_degparms <- backtransform_odeparms(start_degparms,
      orig$mmkin[[1]]$mkinmod,
      transform_rates = orig$mmkin[[1]]$transform_rates,
      transform_fractions = orig$mmkin[[1]]$transform_fractions)
    degparm_names <- names(start_degparms)

    names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index]))

    selected_parms[, degparm_names_transformed] <-
      t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms,
      orig$mmkin[[1]]$mkinmod,
      transform_rates = orig$mmkin[[1]]$transform_rates,
      transform_fractions = orig$mmkin[[1]]$transform_fractions))
    colnames(selected_parms)[1:length(degparm_names)] <- degparm_names
  }

  start_errparms <- orig$so@model@error.init
  names(start_errparms) <- orig$so@model@name.sigma

  start_omegaparms <- orig$so@model@omega.init

  start_parms <- c(start_degparms, start_errparms)

  scale <- match.arg(scale)
  parm_scale <- switch(scale,
    best = selected_parms[which.best(object[selected]), ],
    median = apply(selected_parms, 2, median)
  )

  # Boxplots of all scaled parameters
  selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale))
  boxplot(selected_scaled_parms, log = "y", main = main, ,
    ylab = "Normalised parameters", ...)

  # Show starting parameters
  start_scaled_parms <- rep(NA_real_, length(orig_parms))
  names(start_scaled_parms) <- names(orig_parms)
  start_scaled_parms[names(start_parms)] <-
    start_parms / parm_scale[names(start_parms)]
  points(start_scaled_parms, col = 3, cex = 3)

  # Show parameters of original run
  orig_scaled_parms <- orig_parms / parm_scale
  points(orig_scaled_parms, col = 2, cex = 2)

  abline(h = 1, lty = 2)

  legend(lpos, inset = c(0.05, 0.05), bty = "n",
    pch = 1, col = 3:1, lty = c(NA, NA, 1),
    legend = c(
      "Original start",
      "Original results",
      "Multistart runs"))
}

utils::globalVariables("ds")

#' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
#'
#' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin]
#' @param i A numeric index to select datasets for which to plot the individual predictions,
#' in case plots get too large
#' @inheritParams plot.mkinfit
#' @param standardized Should the residuals be standardized? Only takes effect if
#' `resplot = "time"`.
#' @param pop_curves Per default, one population curve is drawn in case
#' population parameters are fitted by the model, e.g. for saem objects.
#' In case there is a covariate model, the behaviour depends on the value
#' of 'covariates'
#' @param covariates Data frame with covariate values for all variables in
#' any covariate models in the object. If given, it overrides 'covariate_quantiles'.
#' Each line in the data frame will result in a line drawn for the population.
#' Rownames are used in the legend to label the lines.
#' @param covariate_quantiles This argument only has an effect if the fitted
#' object has covariate models. If so, the default is to show three population
#' curves, for the 5th percentile, the 50th percentile and the 95th percentile
#' of the covariate values used for fitting the model.
#' @note Covariate models are currently only supported for saem.mmkin objects.
#' @param pred_over Named list of alternative predictions as obtained
#' from [mkinpredict] with a compatible [mkinmod].
#' @param test_log_parms Passed to [mean_degparms] in the case of an
#' [mixed.mmkin] object
#' @param conf.level Passed to [mean_degparms] in the case of an
#' [mixed.mmkin] object
#' @param default_log_parms Passed to [mean_degparms] in the case of an
#' [mixed.mmkin] object
#' @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 resplot Should the residuals plotted against time or against
#' predicted values?
#' @param col_ds Colors used for plotting the observed data and the
#' corresponding model prediction lines for the different datasets.
#' @param pch_ds Symbols to be used for plotting the data.
#' @param lty_ds Line types to be used for the model predictions.
#' @importFrom stats coefficients
#' @return The function is called for its side effect.
#' @author Johannes Ranke
#' @examples
#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
#'  function(x) x$data[c("name", "time", "value")])
#' names(ds) <- paste0("ds ", 6:10)
#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
#'   A1 = mkinsub("SFO"), quiet = TRUE)
#' \dontrun{
#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
#' plot(f[, 3:4], standardized = TRUE)
#'
#' # For this fit we need to increase pnlsMaxiter, and we increase the
#' # tolerance in order to speed up the fit for this example evaluation
#' # It still takes 20 seconds to run
#' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
#' plot(f_nlme)
#'
#' f_saem <- saem(f, transformations = "saemix")
#' plot(f_saem)
#'
#' f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs")
#' f_nlmix <- nlmix(f_obs)
#' plot(f_nlmix)
#'
#' # We can overlay the two variants if we generate predictions
#' pred_nlme <- mkinpredict(dfop_sfo,
#'   f_nlme$bparms.optim[-1],
#'   c(parent = f_nlme$bparms.optim[[1]], A1 = 0),
#'   seq(0, 180, by = 0.2))
#' plot(f_saem, pred_over = list(nlme = pred_nlme))
#' }
#' @export
plot.mixed.mmkin <- function(x,
  i = 1:ncol(x$mmkin),
  obs_vars = names(x$mkinmod$map),
  standardized = TRUE,
  covariates = NULL,
  covariate_quantiles = c(0.5, 0.05, 0.95),
  xlab = "Time",
  xlim = range(x$data$time),
  resplot = c("predicted", "time"),
  pop_curves = "auto",
  pred_over = NULL,
  test_log_parms = FALSE,
  conf.level = 0.6,
  default_log_parms = NA,
  ymax = "auto", maxabs = "auto",
  ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
  nrow.legend = ceiling((length(i) + 1) / ncol.legend),
  rel.height.legend = 0.02 + 0.07 * nrow.legend,
  rel.height.bottom = 1.1,
  pch_ds = 1:length(i),
  col_ds = pch_ds + 1,
  lty_ds = col_ds,
  frame = TRUE, ...
)
{
  # Prepare parameters and data
  fit_1 <- x$mmkin[[1]]
  ds_names <- colnames(x$mmkin)

  backtransform = TRUE

  if (identical(class(x), "mixed.mmkin")) {
    if (identical(pop_curves, "auto")) {
      pop_curves <- FALSE
    } else {
      pop_curves <- TRUE
    }
    if (pop_curves) {
      degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms,
        conf.level = conf.level, default_log_parms = default_log_parms)
    }

    degparms_tmp <- parms(x$mmkin, transformed = TRUE)
    degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ]))
    residual_type = ifelse(standardized, "standardized", "residual")
    residuals <- x$data[[residual_type]]
  }

  if (inherits(x, "nlme.mmkin")) {
    if (identical(pop_curves, "auto")) {
      pop_curves <- TRUE
    } else {
      pop_curves <- FALSE
    }
    degparms_i <- coefficients(x)
    degparms_pop <- nlme::fixef(x)
    residuals <- residuals(x,
      type = ifelse(standardized, "pearson", "response"))
  }

  if (inherits(x, "saem.mmkin")) {
    if (x$transformations == "saemix") backtransform = FALSE
    psi <- saemix::psi(x$so)
    rownames(psi) <- x$saemix_ds_order
    degparms_i <- psi[ds_names, ]
    degparms_i_names <- colnames(degparms_i)
    residual_type = ifelse(standardized, "standardized", "residual")
    residuals <- x$data[[residual_type]]

    if (identical(pop_curves, "auto")) {
      if (length(x$covariate_models) == 0) {
        degparms_pop <- x$so@results@fixed.effects
        names(degparms_pop) <- degparms_i_names
        pop_curves <- TRUE
      } else {
        if (is.null(covariates)) {
          covariates = as.data.frame(
            apply(x$covariates, 2, quantile,
             covariate_quantiles, simplify = FALSE))
          rownames(covariates) <- paste(
            ifelse(length(x$covariate_models) == 1,
              "Covariate", "Covariates"),
              rownames(covariates))
        }
        degparms_pop <- parms(x, covariates = covariates)
        pop_curves <- TRUE
      }
    } else {
      pop_curves <- FALSE
    }
  }

  if (pop_curves) {
    # Make sure degparms_pop is a matrix, columns corresponding to population curve(s)
    if (is.null(dim(degparms_pop))) {
      degparms_pop <- matrix(degparms_pop, ncol = 1,
        dimnames = list(names(degparms_pop), "Population"))
    }
  }

  degparms_fixed <- fit_1$fixed$value
  names(degparms_fixed) <- rownames(fit_1$fixed)
  degparms_all <- cbind(as.matrix(degparms_i),
    matrix(rep(degparms_fixed, nrow(degparms_i)),
      ncol = length(degparms_fixed),
      nrow = nrow(degparms_i), byrow = TRUE))
  degparms_all_names <- c(names(degparms_i), names(degparms_fixed))
  colnames(degparms_all) <- degparms_all_names

  odeini_names <- grep("_0$", degparms_all_names, value = TRUE)
  odeparms_names <- setdiff(degparms_all_names, odeini_names)

  observed <- cbind(x$data[c("ds", "name", "time", "value")],
    residual = residuals)

  solution_type = fit_1$solution_type

  outtimes <- sort(unique(c(x$data$time,
    seq(xlim[1], xlim[2], length.out = 50))))

  pred_list <- lapply(i, function(ds_i)   {
    odeparms_trans <- degparms_all[ds_i, odeparms_names]
    names(odeparms_trans) <- odeparms_names # needed if only one odeparm
    if (backtransform) {
      odeparms <- backtransform_odeparms(odeparms_trans,
        x$mkinmod,
        transform_rates = fit_1$transform_rates,
        transform_fractions = fit_1$transform_fractions)
    } else {
      odeparms <- odeparms_trans
    }

    odeini <- degparms_all[ds_i, odeini_names]
    names(odeini) <- gsub("_0", "", odeini_names)

    out <- mkinpredict(x$mkinmod, odeparms, odeini,
      outtimes, solution_type = solution_type,
      atol = fit_1$atol, rtol = fit_1$rtol)
  })
  names(pred_list) <- ds_names[i]
  pred_ds <- vctrs::vec_rbind(!!!pred_list, .names_to = "ds")

  if (pop_curves) {
    pred_list_pop <- lapply(1:ncol(degparms_pop), function(cov_i)   {
      degparms_all_pop_i <- c(degparms_pop[, cov_i], degparms_fixed)
      odeparms_pop_trans_i <- degparms_all_pop_i[odeparms_names]
      names(odeparms_pop_trans_i) <- odeparms_names # needed if only one odeparm
      if (backtransform) {
        odeparms_pop_i <- backtransform_odeparms(odeparms_pop_trans_i,
          x$mkinmod,
          transform_rates = fit_1$transform_rates,
          transform_fractions = fit_1$transform_fractions)
      } else {
        odeparms_pop_i <- odeparms_pop_trans_i
      }

      odeini <- degparms_all_pop_i[odeini_names]
      names(odeini) <- gsub("_0", "", odeini_names)

      out <- mkinpredict(x$mkinmod, odeparms_pop_i, odeini,
        outtimes, solution_type = solution_type,
        atol = fit_1$atol, rtol = fit_1$rtol)
    })
    names(pred_list_pop) <- colnames(degparms_pop)

  } else {
    pred_list_pop <- NULL
  }

  # Start of graphical section
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar, no.readonly = TRUE))

  n_plot_rows = length(obs_vars)
  n_plots = n_plot_rows * 2

  # Set relative plot heights, so the first plot row is the norm
  rel.heights <- if (n_plot_rows > 1) {
    c(rel.height.legend, c(rep(1, n_plot_rows - 1), rel.height.bottom))
  } else {
    c(rel.height.legend, 1)
  }

  layout_matrix = matrix(c(1, 1, 2:(n_plots + 1)),
    n_plot_rows + 1, 2, byrow = TRUE)
  layout(layout_matrix, heights = rel.heights)

  par(mar = c(0.1, 2.1, 0.1, 2.1))

  # Empty plot with legend
  if (!is.null(pred_over)) lty_over <- seq(2, length.out = length(pred_over))
  else lty_over <- NULL
  if (pop_curves) {
    if (is.null(covariates)) {
      lty_pop <- 1
      names(lty_pop) <- "Population"
    } else {
      lty_pop <- 1:nrow(covariates)
      names(lty_pop) <- rownames(covariates)
    }
  } else {
    lty_pop <- NULL
  }
  n_pop_over <- length(lty_pop) + length(lty_over)

  plot(0, type = "n", axes = FALSE, ann = FALSE)
  legend("center", bty = "n", ncol = ncol.legend,
    legend = c(names(lty_pop), names(pred_over), ds_names[i]),
    lty = c(lty_pop, lty_over, lty_ds),
    lwd = c(rep(2, n_pop_over), rep(1, length(i))),
    col = c(rep(1, n_pop_over), col_ds),
    pch = c(rep(NA, n_pop_over), pch_ds))

  resplot <- match.arg(resplot)

  # Loop plot rows
  for (plot_row in 1:n_plot_rows) {

    obs_var <- obs_vars[plot_row]
    observed_row <- subset(observed, name == obs_var)

    # Set ylim to sensible default, or use ymax
    if (identical(ymax, "auto")) {
      ylim_row = c(0,
        max(c(observed_row$value, pred_ds[[obs_var]]), na.rm = TRUE))
    } else {
      ylim_row = c(0, ymax[plot_row])
    }

    # Margins for bottom row of plots when we have more than one row
    # This is the only row that needs to show the x axis legend
    if (plot_row == n_plot_rows) {
      par(mar = c(5.1, 4.1, 1.1, 2.1))
    } else {
      par(mar = c(3.0, 4.1, 1.1, 2.1))
    }

    plot(0, type = "n",
      xlim = xlim, ylim = ylim_row,
      xlab = xlab, ylab = paste("Residues", obs_var), frame = frame)

    if (!is.null(pred_over)) {
      for (i_over in seq_along(pred_over)) {
        pred_frame <- as.data.frame(pred_over[[i_over]])
        lines(pred_frame$time, pred_frame[[obs_var]],
          lwd = 2, lty = lty_over[i_over])
      }
    }

    for (ds_i in seq_along(i)) {
      points(subset(observed_row, ds == ds_names[ds_i], c("time", "value")),
        col = col_ds[ds_i], pch = pch_ds[ds_i])
      lines(subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)),
        col = col_ds[ds_i], lty = lty_ds[ds_i])
    }

    if (pop_curves) {
      for (cov_i in seq_along(pred_list_pop)) {
        cov_name <- names(pred_list_pop)[cov_i]
        lines(
          pred_list_pop[[cov_i]][, "time"],
          pred_list_pop[[cov_i]][, obs_var],
          type = "l", lwd = 2, lty = lty_pop[cov_i])
      }
    }

    if (identical(maxabs, "auto")) {
      maxabs = max(abs(observed_row$residual), na.rm = TRUE)
    }

    if (identical(resplot, "time")) {
      plot(0, type = "n", xlim = xlim, xlab = "Time",
        ylim = c(-1.2 * maxabs, 1.2 * maxabs),
        ylab = if (standardized) "Standardized residual" else "Residual",
        frame = frame)

      abline(h = 0, lty = 2)

      for (ds_i in seq_along(i)) {
        points(subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
          col = col_ds[ds_i], pch = pch_ds[ds_i])
      }
    }

    if (identical(resplot, "predicted")) {
      plot(0, type = "n",
        xlim = c(0, max(pred_ds[[obs_var]])),
        xlab = "Predicted",
        ylim = c(-1.2 * maxabs, 1.2 * maxabs),
        ylab = if (standardized) "Standardized residual" else "Residual",
        frame = frame)

      abline(h = 0, lty = 2)

      for (ds_i in seq_along(i)) {
        observed_row_ds <- merge(
          subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
          subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)))
        points(observed_row_ds[c(3, 2)],
          col = col_ds[ds_i], pch = pch_ds[ds_i])
      }
    }
  }
}

#' Summary method for class "mkinfit"
#'
#' Lists model equations, initial parameter values, optimised parameters with
#' some uncertainty statistics, the chi2 error levels calculated according to
#' FOCUS guidance (2006) as defined therein, formation fractions, DT50 values
#' and optionally the data, consisting of observed, predicted and residual
#' values.
#'
#' @param object an object of class [mkinfit].
#' @param x an object of class \code{summary.mkinfit}.
#' @param data logical, indicating whether the data should be included in the
#' summary.
#' @param distimes logical, indicating whether DT50 and DT90 values should be
#' included.
#' @param alpha error level for confidence interval estimation from t
#' distribution
#' @param digits Number of digits to use for printing
#' @param \dots optional arguments passed to methods like \code{print}.
#' @importFrom stats qt pt cov2cor
#' @return The summary function returns a list with components, among others
#'   \item{version, Rversion}{The 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 model}
#'   \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
#'   \item{bpar}{Optimised and backtransformed
#'     parameters}
#'   \item{data}{The data (see Description above).}
#'   \item{start}{The starting values and bounds, if applicable, for optimised
#'     parameters.}
#'   \item{fixed}{The values of fixed parameters.}
#'   \item{errmin }{The chi2 error levels for
#'     each observed variable.}
#'   \item{bparms.ode}{All backtransformed ODE
#'     parameters, for use as starting parameters for related models.}
#'   \item{errparms}{Error model parameters.}
#'   \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 and fractional eigenvector component
#'      g of SFORB systems in the model.}
#'   The print method is called for its side effect, i.e. printing the summary.
#' @author Johannes Ranke
#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence
#'   and Degradation Kinetics from Environmental Fate Studies on Pesticides in
#'   EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
#'   EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
#'   \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
#' @examples
#'
#'   summary(mkinfit("SFO", FOCUS_2006_A, quiet = TRUE))
#'
#' @export
summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
  param  <- object$par
  pnames <- names(param)
  bpnames <- names(object$bparms.optim)
  epnames <- names(object$errparms)
  p      <- length(param)
  mod_vars <- names(object$mkinmod$diffs)
  covar  <- try(solve(object$hessian), silent = TRUE)
  covar_notrans  <- try(solve(object$hessian_notrans), silent = TRUE)
  rdf <- object$df.residual

  if (!is.numeric(covar) | is.na(covar[1])) {
    covar <- NULL
    se <- lci <- uci <- rep(NA, p)
  } else {
    rownames(covar) <- colnames(covar) <- pnames
    se     <- sqrt(diag(covar))
    lci    <- param + qt(alpha/2, rdf) * se
    uci    <- param + qt(1-alpha/2, rdf) * se
  }

  beparms.optim <- c(object$bparms.optim, object$par[epnames])
  if (!is.numeric(covar_notrans) | is.na(covar_notrans[1])) {
    covar_notrans <- NULL
    se_notrans <- tval <- pval <- rep(NA, p)
  } else {
    rownames(covar_notrans) <- colnames(covar_notrans) <- c(bpnames, epnames)
    se_notrans <- sqrt(diag(covar_notrans))
    tval  <- beparms.optim / se_notrans
    pval  <- pt(abs(tval), rdf, lower.tail = FALSE)
  }

  names(se) <- pnames

  param <- cbind(param, se, lci, uci)
  dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))

  bparam <- cbind(Estimate = beparms.optim, se_notrans,
                  "t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA)

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

  for (pname in pnames) {
    if (!pname %in% f_names_skip) {
      par.lower <- param[pname, "Lower"]
      par.upper <- param[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)
      bparam[names(bpl), "Lower"] <- bpl
      bparam[names(bpu), "Upper"] <- bpu
    }
  }
  bparam[epnames, c("Lower", "Upper")] <- param[epnames, c("Lower", "Upper")]

  ans <- list(
    version = as.character(utils::packageVersion("mkin")),
    Rversion = paste(R.version$major, R.version$minor, sep="."),
    date.fit = object$date,
    date.summary = date(),
    solution_type = object$solution_type,
    warnings = object$summary_warnings,
    use_of_ff = object$mkinmod$use_of_ff,
    error_model_algorithm = object$error_model_algorithm,
    df = c(p, rdf),
    covar = covar,
    covar_notrans = covar_notrans,
    err_mod = object$err_mod,
    niter = object$iterations,
    calls = object$calls,
    time = object$time,
    par = param,
    bpar = bparam)

  if (!is.null(object$version)) {
    ans$fit_version <- object$version
    ans$fit_Rversion <- object$Rversion
    if (ans$fit_version >= "0.9.49.6") {
      ans$AIC = AIC(object)
      ans$BIC = BIC(object)
      ans$logLik = logLik(object)
    }
  }

  ans$diffs <- object$mkinmod$diffs
  if(data) ans$data <- object$data
  ans$start <- object$start
  ans$start_transformed <- object$start_transformed

  ans$fixed <- object$fixed

  ans$errmin <- mkinerrmin(object, alpha = 0.05)

  if (object$calls > 0) {
    if (!is.null(ans$covar)){
      Corr <- cov2cor(ans$covar)
      rownames(Corr) <- colnames(Corr) <- rownames(ans$par)
      ans$Corr <- Corr
    } else {
      warning("Could not calculate correlation; no covariance matrix")
    }
  }

  ans$bparms.ode <- object$bparms.ode
  ans$shapiro.p <- object$shapiro.p
  ep <- endpoints(object)
  if (length(ep$ff) != 0)
    ans$ff <- ep$ff
  if (distimes) ans$distimes <- ep$distimes
  if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
  if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message
  class(ans) <- "summary.mkinfit"
  return(ans)
}

#' @rdname summary.mkinfit
#' @export
print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
  if (is.null(x$fit_version)) {
    cat("mkin version:   ", x$version, "\n")
    cat("R version:      ", x$Rversion, "\n")
  } else {
    cat("mkin version used for fitting:   ", x$fit_version, "\n")
    cat("R version used for fitting:      ", x$fit_Rversion, "\n")
  }

  cat("Date of fit:    ", x$date.fit, "\n")
  cat("Date of summary:", x$date.summary, "\n")


1		#' Plot parameter variability of multistart objects
2		#'
3		#' Produces a boxplot with all parameters from the multiple runs, scaled
4		#' either by the parameters of the run with the highest likelihood,
5		#' or by their medians as proposed in the paper by Duchesne et al. (2021).
6		#'
7		#' Starting values of degradation model parameters and error model parameters
8		#' are shown as green circles. The results obtained in the original run
9		#' are shown as red circles.
10		#'
11		#' @param object The [multistart] object
12		#' @param llmin The minimum likelihood of objects to be shown
13		#' @param llquant Fractional value for selecting only the fits with higher
14		#' likelihoods. Overrides 'llmin'.
15		#' @param scale By default, scale parameters using the best
16		#' available fit.
17		#' If 'median', parameters are scaled using the median parameters from all fits.
18		#' @param main Title of the plot
19		#' @param lpos Positioning of the legend.
20		#' @param \dots Passed to [boxplot]
21		#' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical
22		#' identifiability in the frame of nonlinear mixed effects models: the example
23		#' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478.
24		#' doi: 10.1186/s12859-021-04373-4.
25		#' @seealso [multistart]
26		#' @importFrom stats median quantile
27		#' @export
28		parplot <- function(object, ...) {
29	176x	UseMethod("parplot")
30		}
31
32		#' @rdname parplot
33		#' @export
34		parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
35		scale = c("best", "median"),
36		lpos = "bottomleft", main = "", ...)
37		{
38	176x	oldpar <- par(no.readonly = TRUE)
39	176x	on.exit(par(oldpar, no.readonly = TRUE))
40
41	176x	orig <- attr(object, "orig")
42	176x	orig_parms <- parms(orig)
43	176x	start_degparms <- orig$mean_dp_start
44	176x	all_parms <- parms(object, exclude_failed = FALSE)
45
46	176x	if (inherits(object, "multistart.saem.mmkin")) {
47	176x	llfunc <- function(object) {
48	!	if (inherits(object$so, "try-error")) return(NA)
49	1408x	else return(logLik(object$so))
50		}
51		} else {
52	!	stop("parplot is only implemented for multistart.saem.mmkin objects")
53		}
54	176x	ll <- sapply(object, llfunc)
55	176x	if (!is.na(llquant[1])) {
56	!	if (llmin != -Inf) warning("Overriding 'llmin' because 'llquant' was specified")
57	88x	llmin <- quantile(ll, 1 - llquant)
58		}
59	176x	selected <- which(ll > llmin)
60	176x	selected_parms <- all_parms[selected, ]
61
62	176x	par(las = 1)
63	176x	if (orig$transformations == "mkin") {
64	88x	degparm_names_transformed <- names(start_degparms)
65	88x	degparm_index <- which(names(orig_parms) %in% degparm_names_transformed)
66	88x	orig_parms[degparm_names_transformed] <- backtransform_odeparms(
67	88x	orig_parms[degparm_names_transformed],
68	88x	orig$mmkin[[1]]$mkinmod,
69	88x	transform_rates = orig$mmkin[[1]]$transform_rates,
70	88x	transform_fractions = orig$mmkin[[1]]$transform_fractions)
71	88x	start_degparms <- backtransform_odeparms(start_degparms,
72	88x	orig$mmkin[[1]]$mkinmod,
73	88x	transform_rates = orig$mmkin[[1]]$transform_rates,
74	88x	transform_fractions = orig$mmkin[[1]]$transform_fractions)
75	88x	degparm_names <- names(start_degparms)
76
77	88x	names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index]))
78
79	88x	selected_parms[, degparm_names_transformed] <-
80	88x	t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms,
81	88x	orig$mmkin[[1]]$mkinmod,
82	88x	transform_rates = orig$mmkin[[1]]$transform_rates,
83	88x	transform_fractions = orig$mmkin[[1]]$transform_fractions))
84	88x	colnames(selected_parms)[1:length(degparm_names)] <- degparm_names
85		}
86
87	176x	start_errparms <- orig$so@model@error.init
88	176x	names(start_errparms) <- orig$so@model@name.sigma
89
90	176x	start_omegaparms <- orig$so@model@omega.init
91
92	176x	start_parms <- c(start_degparms, start_errparms)
93
94	176x	scale <- match.arg(scale)
95	176x	parm_scale <- switch(scale,
96	176x	best = selected_parms[which.best(object[selected]), ],
97	176x	median = apply(selected_parms, 2, median)
98		)
99
100		# Boxplots of all scaled parameters
101	176x	selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale))
102	176x	boxplot(selected_scaled_parms, log = "y", main = main, ,
103	176x	ylab = "Normalised parameters", ...)
104
105		# Show starting parameters
106	176x	start_scaled_parms <- rep(NA_real_, length(orig_parms))
107	176x	names(start_scaled_parms) <- names(orig_parms)
108	176x	start_scaled_parms[names(start_parms)] <-
109	176x	start_parms / parm_scale[names(start_parms)]
110	176x	points(start_scaled_parms, col = 3, cex = 3)
111
112		# Show parameters of original run
113	176x	orig_scaled_parms <- orig_parms / parm_scale
114	176x	points(orig_scaled_parms, col = 2, cex = 2)
115
116	176x	abline(h = 1, lty = 2)
117
118	176x	legend(lpos, inset = c(0.05, 0.05), bty = "n",
119	176x	pch = 1, col = 3:1, lty = c(NA, NA, 1),
120	176x	legend = c(
121	176x	"Original start",
122	176x	"Original results",
123	176x	"Multistart runs"))
124		}

1		#' Summary method for class "mkinfit"
2		#'
3		#' Lists model equations, initial parameter values, optimised parameters with
4		#' some uncertainty statistics, the chi2 error levels calculated according to
5		#' FOCUS guidance (2006) as defined therein, formation fractions, DT50 values
6		#' and optionally the data, consisting of observed, predicted and residual
7		#' values.
8		#'
9		#' @param object an object of class [mkinfit].
10		#' @param x an object of class \code{summary.mkinfit}.
11		#' @param data logical, indicating whether the data should be included in the
12		#' summary.
13		#' @param distimes logical, indicating whether DT50 and DT90 values should be
14		#' included.
15		#' @param alpha error level for confidence interval estimation from t
16		#' distribution
17		#' @param digits Number of digits to use for printing
18		#' @param \dots optional arguments passed to methods like \code{print}.
19		#' @importFrom stats qt pt cov2cor
20		#' @return The summary function returns a list with components, among others
21		#' \item{version, Rversion}{The mkin and R versions used}
22		#' \item{date.fit, date.summary}{The dates where the fit and the summary were
23		#' produced}
24		#' \item{diffs}{The differential equations used in the model}
25		#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
26		#' \item{bpar}{Optimised and backtransformed
27		#' parameters}
28		#' \item{data}{The data (see Description above).}
29		#' \item{start}{The starting values and bounds, if applicable, for optimised
30		#' parameters.}
31		#' \item{fixed}{The values of fixed parameters.}
32		#' \item{errmin }{The chi2 error levels for
33		#' each observed variable.}
34		#' \item{bparms.ode}{All backtransformed ODE
35		#' parameters, for use as starting parameters for related models.}
36		#' \item{errparms}{Error model parameters.}
37		#' \item{ff}{The estimated formation fractions derived from the fitted
38		#' model.}
39		#' \item{distimes}{The DT50 and DT90 values for each observed variable.}
40		#' \item{SFORB}{If applicable, eigenvalues and fractional eigenvector component
41		#' g of SFORB systems in the model.}
42		#' The print method is called for its side effect, i.e. printing the summary.
43		#' @author Johannes Ranke
44		#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence
45		#' and Degradation Kinetics from Environmental Fate Studies on Pesticides in
46		#' EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
47		#' EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
48		#' \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
49		#' @examples
50		#'
51		#' summary(mkinfit("SFO", FOCUS_2006_A, quiet = TRUE))
52		#'
53		#' @export
54		summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
55	52158x	param <- object$par
56	52158x	pnames <- names(param)
57	52158x	bpnames <- names(object$bparms.optim)
58	52158x	epnames <- names(object$errparms)
59	52158x	p <- length(param)
60	52158x	mod_vars <- names(object$mkinmod$diffs)
61	52158x	covar <- try(solve(object$hessian), silent = TRUE)
62	52158x	covar_notrans <- try(solve(object$hessian_notrans), silent = TRUE)
63	52158x	rdf <- object$df.residual
64
65	52158x	if (!is.numeric(covar) \| is.na(covar[1])) {
66	!	covar <- NULL
67	!	se <- lci <- uci <- rep(NA, p)
68		} else {
69	52158x	rownames(covar) <- colnames(covar) <- pnames
70	52158x	se <- sqrt(diag(covar))
71	52158x	lci <- param + qt(alpha/2, rdf) * se
72	52158x	uci <- param + qt(1-alpha/2, rdf) * se
73		}
74
75	52158x	beparms.optim <- c(object$bparms.optim, object$par[epnames])
76	52158x	if (!is.numeric(covar_notrans) \| is.na(covar_notrans[1])) {
77	88x	covar_notrans <- NULL
78	88x	se_notrans <- tval <- pval <- rep(NA, p)
79		} else {
80	52070x	rownames(covar_notrans) <- colnames(covar_notrans) <- c(bpnames, epnames)
81	52070x	se_notrans <- sqrt(diag(covar_notrans))
82	52070x	tval <- beparms.optim / se_notrans
83	52070x	pval <- pt(abs(tval), rdf, lower.tail = FALSE)
84		}
85
86	52158x	names(se) <- pnames
87
88	52158x	param <- cbind(param, se, lci, uci)
89	52158x	dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))
90
91	52158x	bparam <- cbind(Estimate = beparms.optim, se_notrans,
92	52158x	"t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA)
93
94		# Transform boundaries of CI for one parameter at a time,
95		# with the exception of sets of formation fractions (single fractions are OK).
96	52158x	f_names_skip <- character(0)
97	52158x	for (box in mod_vars) { # Figure out sets of fractions to skip
98	70671x	f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
99	70671x	n_paths <- length(f_names)
100	1135x	if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
101		}
102
103	52158x	for (pname in pnames) {
104	293621x	if (!pname %in% f_names_skip) {
105	290217x	par.lower <- param[pname, "Lower"]
106	290217x	par.upper <- param[pname, "Upper"]
107	290217x	names(par.lower) <- names(par.upper) <- pname
108	290217x	bpl <- backtransform_odeparms(par.lower, object$mkinmod,
109	290217x	object$transform_rates,
110	290217x	object$transform_fractions)
111	290217x	bpu <- backtransform_odeparms(par.upper, object$mkinmod,
112	290217x	object$transform_rates,
113	290217x	object$transform_fractions)
114	290217x	bparam[names(bpl), "Lower"] <- bpl
115	290217x	bparam[names(bpu), "Upper"] <- bpu
116		}
117		}
118	52158x	bparam[epnames, c("Lower", "Upper")] <- param[epnames, c("Lower", "Upper")]
119
120	52158x	ans <- list(
121	52158x	version = as.character(utils::packageVersion("mkin")),
122	52158x	Rversion = paste(R.version$major, R.version$minor, sep="."),
123	52158x	date.fit = object$date,
124	52158x	date.summary = date(),
125	52158x	solution_type = object$solution_type,
126	52158x	warnings = object$summary_warnings,
127	52158x	use_of_ff = object$mkinmod$use_of_ff,
128	52158x	error_model_algorithm = object$error_model_algorithm,
129	52158x	df = c(p, rdf),
130	52158x	covar = covar,
131	52158x	covar_notrans = covar_notrans,
132	52158x	err_mod = object$err_mod,
133	52158x	niter = object$iterations,
134	52158x	calls = object$calls,
135	52158x	time = object$time,
136	52158x	par = param,
137	52158x	bpar = bparam)
138
139	52158x	if (!is.null(object$version)) {
140	52158x	ans$fit_version <- object$version
141	52158x	ans$fit_Rversion <- object$Rversion
142	52158x	if (ans$fit_version >= "0.9.49.6") {
143	52156x	ans$AIC = AIC(object)
144	52156x	ans$BIC = BIC(object)
145	52156x	ans$logLik = logLik(object)
146		}
147		}
148
149	52158x	ans$diffs <- object$mkinmod$diffs
150	52158x	if(data) ans$data <- object$data
151	52158x	ans$start <- object$start
152	52158x	ans$start_transformed <- object$start_transformed
153
154	52158x	ans$fixed <- object$fixed
155
156	52158x	ans$errmin <- mkinerrmin(object, alpha = 0.05)
157
158	52158x	if (object$calls > 0) {
159	52158x	if (!is.null(ans$covar)){
160	52158x	Corr <- cov2cor(ans$covar)
161	52158x	rownames(Corr) <- colnames(Corr) <- rownames(ans$par)
162	52158x	ans$Corr <- Corr
163		} else {
164	!	warning("Could not calculate correlation; no covariance matrix")
165		}
166		}
167
168	52158x	ans$bparms.ode <- object$bparms.ode
169	52158x	ans$shapiro.p <- object$shapiro.p
170	52158x	ep <- endpoints(object)
171	52158x	if (length(ep$ff) != 0)
172	15612x	ans$ff <- ep$ff
173	52158x	if (distimes) ans$distimes <- ep$distimes
174	2442x	if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
175	43972x	if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message
176	52158x	class(ans) <- "summary.mkinfit"
177	52158x	return(ans)
178		}
179
180		#' @rdname summary.mkinfit
181		#' @export
182		print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
183	4x	if (is.null(x$fit_version)) {
184	!	cat("mkin version: ", x$version, "\n")
185	!	cat("R version: ", x$Rversion, "\n")
186		} else {
187	4x	cat("mkin version used for fitting: ", x$fit_version, "\n")
188	4x	cat("R version used for fitting: ", x$fit_Rversion, "\n")
189		}
190
191	4x	cat("Date of fit: ", x$date.fit, "\n")
192	4x	cat("Date of summary:", x$date.summary, "\n")

1		#' Method to get status information for fit array objects
2		#'
3		#' @param object The object to investigate
4		#' @param x The object to be printed
5		#' @param \dots For potential future extensions
6		#' @return An object with the same dimensions as the fit array
7		#' suitable printing method.
8		#' @export
9		status <- function(object, ...)
10		{
11	589x	UseMethod("status", object)
12		}
13
14		#' @rdname status
15		#' @export
16		#' @examples
17		#' \dontrun{
18		#' fits <- mmkin(
19		#' c("SFO", "FOMC"),
20		#' list("FOCUS A" = FOCUS_2006_A,
21		#' "FOCUS B" = FOCUS_2006_C),
22		#' quiet = TRUE)
23		#' status(fits)
24		#' }
25		status.mmkin <- function(object, ...) {
26	376x	all_summary_warnings <- character()
27	376x	sww <- 0 # Counter for Shapiro-Wilks warnings
28
29	376x	result <- lapply(object,
30	376x	function(fit) {
31	!	if (inherits(fit, "try-error")) return("E")
32	4391x	sw <- fit$summary_warnings
33	4391x	swn <- names(sw)
34	4391x	if (length(sw) > 0) {
35	!	if (any(grepl("S", swn))) {
36	!	sww <<- sww + 1
37	!	swn <- gsub("S", paste0("S", sww), swn)
38		}
39	!	warnstring <- paste(swn, collapse = ", ")
40	!	names(sw) <- swn
41	!	all_summary_warnings <<- c(all_summary_warnings, sw)
42	!	return(warnstring)
43		} else {
44	4391x	return("OK")
45		}
46		})
47	376x	result <- unlist(result)
48	376x	dim(result) <- dim(object)
49	376x	dimnames(result) <- dimnames(object)
50
51	376x	u_swn <- unique(names(all_summary_warnings))
52	376x	attr(result, "unique_warnings") <- all_summary_warnings[u_swn]
53	376x	class(result) <- "status.mmkin"
54	376x	return(result)
55		}
56
57		#' @rdname status
58		#' @export
59		print.status.mmkin <- function(x, ...) {
60	376x	u_w <- attr(x, "unique_warnings")
61	376x	attr(x, "unique_warnings") <- NULL
62	376x	class(x) <- NULL
63	376x	print(x, quote = FALSE)
64	376x	cat("\n")
65	376x	for (i in seq_along(u_w)) {
66	!	cat(names(u_w)[i], ": ", u_w[i], "\n", sep = "")
67		}
68	376x	if (any(x == "OK")) cat("OK: No warnings\n")
69	!	if (any(x == "E")) cat("E: Error\n")
70		}
71
72		#' @rdname status
73		#' @export
74		status.mhmkin <- function(object, ...) {
75	125x	if (inherits(object[[1]], "saem.mmkin")) {
76	125x	test_func <- function(fit) {
77	500x	if (inherits(fit, "try-error")) {
78	!	return("E")
79		} else {
80	500x	if (inherits(fit$so, "try-error")) {
81	!	return("E")
82		} else {
83	500x	if (!is.null(fit$FIM_failed)) {
84	!	return_values <- c("fixed effects" = "Fth",
85	!	"random effects and error model parameters" = "FO")
86	!	return(paste(return_values[fit$FIM_failed], collapse = ", "))
87		} else {
88	500x	return("OK")
89		}
90		}
91		}
92		}
93		} else {
94	!	stop("Only mhmkin objects containing saem.mmkin objects currently supported")
95		}
96	125x	result <- lapply(object, test_func)
97	125x	result <- unlist(result)
98	125x	dim(result) <- dim(object)
99	125x	dimnames(result) <- dimnames(object)
100
101	125x	class(result) <- "status.mhmkin"
102	125x	return(result)
103		}
104
105		#' @rdname status
106		#' @export
107		print.status.mhmkin <- function(x, ...) {
108	125x	class(x) <- NULL
109	125x	print(x, quote = FALSE)
110	125x	cat("\n")
111	125x	if (any(x == "OK")) cat("OK: Fit terminated successfully\n")
112	!	if (any(x == "Fth")) cat("Fth: Could not invert FIM for fixed effects\n")
113	!	if (any(x == "FO")) cat("FO: Could not invert FIM for random effects and error model parameters\n")
114	!	if (any(x == "Fth, FO")) cat("Fth, FO: Could not invert FIM for fixed effects, nor for random effects and error model parameters\n")
115	!	if (any(x == "E")) cat("E: Error\n")
116		}
117

1		#' Function to calculate endpoints for further use from kinetic models fitted
2		#' with mkinfit
3		#'
4		#' This function calculates DT50 and DT90 values as well as formation fractions
5		#' from kinetic models fitted with mkinfit. If the SFORB model was specified
6		#' for one of the parents or metabolites, the Eigenvalues are returned. These
7		#' are equivalent to the rate constants of the DFOP model, but with the
8		#' advantage that the SFORB model can also be used for metabolites.
9		#'
10		#' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
11		#' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
12		#'
13		#' @param fit An object of class [mkinfit], [nlme.mmkin] or [saem.mmkin], or
14		#' another object that has list components mkinmod containing an [mkinmod]
15		#' degradation model, and two numeric vectors, bparms.optim and bparms.fixed,
16		#' that contain parameter values for that model.
17		#' @param covariates Numeric vector with covariate values for all variables in
18		#' any covariate models in the object. If given, it overrides 'covariate_quantile'.
19		#' @param covariate_quantile This argument only has an effect if the fitted
20		#' object has covariate models. If so, the default is to show endpoints
21		#' for the median of the covariate values (50th percentile).
22		#' @importFrom stats optimize
23		#' @return A list with a matrix of dissipation times named distimes, and, if
24		#' applicable, a vector of formation fractions named ff and, if the SFORB model
25		#' was in use, a vector of eigenvalues of these SFORB models, equivalent to
26		#' DFOP rate constants
27		#' @note The function is used internally by [summary.mkinfit],
28		#' [summary.nlme.mmkin] and [summary.saem.mmkin].
29		#' @author Johannes Ranke
30		#' @examples
31		#'
32		#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
33		#' endpoints(fit)
34		#' \dontrun{
35		#' fit_2 <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE)
36		#' endpoints(fit_2)
37		#' fit_3 <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE)
38		#' endpoints(fit_3)
39		#' }
40		#'
41		#' @export
42		endpoints <- function(fit, covariates = NULL, covariate_quantile = 0.5) {
43	56208x	mkinmod <- fit$mkinmod
44	56208x	obs_vars <- names(mkinmod$spec)
45
46	56208x	if (!is.null(fit$covariate_models)) {
47	110x	if (is.null(covariates)) {
48	110x	covariates = as.data.frame(
49	110x	apply(fit$covariates, 2, quantile,
50	110x	covariate_quantile, simplify = FALSE))
51		} else {
52	!	covariate_m <- matrix(covariates, byrow = TRUE)
53	!	colnames(covariate_m) <- names(covariates)
54	!	rownames(covariate_m) <- "User"
55	!	covariates <- as.data.frame(covariate_m)
56		}
57	110x	degparms_trans <- parms(fit, covariates = covariates)[, 1]
58	110x	if (inherits(fit, "saem.mmkin") & (fit$transformations == "saemix")) {
59	!	degparms <- degparms_trans
60		} else {
61	110x	degparms <- backtransform_odeparms(degparms_trans,
62	110x	fit$mkinmod,
63	110x	transform_rates = fit$transform_rates,
64	110x	transform_fractions = fit$transform_fractions)
65		}
66		} else {
67	56098x	degparms <- c(fit$bparms.optim, fit$bparms.fixed)
68		}
69
70		# Set up object to return
71	56208x	ep <- list()
72	56208x	ep$covariates <- covariates
73	56208x	ep$ff <- vector()
74	56208x	ep$SFORB <- vector()
75	56208x	ep$distimes <- data.frame(
76	56208x	DT50 = rep(NA, length(obs_vars)),
77	56208x	DT90 = rep(NA, length(obs_vars)),
78	56208x	row.names = obs_vars)
79
80	56208x	for (obs_var in obs_vars) {
81	73858x	type = names(mkinmod$map[[obs_var]])[1]
82
83		# Get formation fractions if directly fitted, and calculate remaining fraction to sink
84	73858x	f_names = grep(paste("^f", obs_var, sep = "_"), names(degparms), value=TRUE)
85	73858x	if (length(f_names) > 0) {
86	15068x	f_values = degparms[f_names]
87	15068x	f_to_sink = 1 - sum(f_values)
88	15068x	names(f_to_sink) = ifelse(type == "SFORB",
89	15068x	paste(obs_var, "free", "sink", sep = "_"),
90	15068x	paste(obs_var, "sink", sep = "_"))
91	15068x	for (f_name in f_names) {
92	17338x	ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
93		}
94	15068x	ep$ff = append(ep$ff, f_to_sink)
95		}
96
97		# Get the rest
98	73858x	if (type == "SFO") {
99	40900x	k_names = grep(paste("^k", obs_var, sep="_"), names(degparms), value=TRUE)
100	40900x	k_tot = sum(degparms[k_names])
101	40900x	DT50 = log(2)/k_tot
102	40900x	DT90 = log(10)/k_tot
103	40900x	if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) {
104	622x	for (k_name in k_names)
105		{
106	932x	ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
107		}
108		}
109		}
110	73858x	if (type == "FOMC") {
111	1790x	alpha = degparms["alpha"]
112	1790x	beta = degparms["beta"]
113	1790x	DT50 = beta * (2^(1/alpha) - 1)
114	1790x	DT90 = beta * (10^(1/alpha) - 1)
115	1790x	DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
116	1790x	ep$distimes[obs_var, c("DT50back")] = DT50_back
117		}
118	73858x	if (type == "IORE") {
119	352x	k_names = grep(paste("^k__iore", obs_var, sep="_"), names(degparms), value=TRUE)
120	352x	k_tot = sum(degparms[k_names])
121		# From the NAFTA kinetics guidance, p. 5
122	352x	n = degparms[paste("N", obs_var, sep = "_")]
123	352x	k = k_tot
124		# Use the initial concentration of the parent compound
125	352x	source_name = mkinmod$map[[1]][[1]]
126	352x	c0 = degparms[paste(source_name, "0", sep = "_")]
127	352x	alpha = 1 / (n - 1)
128	352x	beta = (c0^(1 - n))/(k * (n - 1))
129	352x	DT50 = beta * (2^(1/alpha) - 1)
130	352x	DT90 = beta * (10^(1/alpha) - 1)
131	352x	DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
132	352x	ep$distimes[obs_var, c("DT50back")] = DT50_back
133	352x	if (mkinmod$use_of_ff == "min") {
134	!	for (k_name in k_names)
135		{
136	!	ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
137		}
138		}
139		}
140	73858x	if (type == "DFOP") {
141	27729x	k1 = degparms["k1"]
142	27729x	k2 = degparms["k2"]
143	27729x	g = degparms["g"]
144	27729x	f <- function(log_t, x) {
145	684705x	t <- exp(log_t)
146	684705x	fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
147	684705x	(fraction - (1 - x/100))^2
148		}
149	27729x	DT50_k1 = log(2)/k1
150	27729x	DT50_k2 = log(2)/k2
151	27729x	DT90_k1 = log(10)/k1
152	27729x	DT90_k2 = log(10)/k2
153
154	27729x	DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum),
155	27729x	silent = TRUE)
156	27729x	DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum),
157	27729x	silent = TRUE)
158	!	if (inherits(DT50, "try-error")) DT50 = NA
159	!	if (inherits(DT90, "try-error")) DT90 = NA
160	27729x	DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
161
162	27729x	ep$distimes[obs_var, c("DT50back")] = DT50_back
163	27729x	ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
164	27729x	ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
165		}
166	73858x	if (type == "HS") {
167	318x	k1 = degparms["k1"]
168	318x	k2 = degparms["k2"]
169	318x	tb = degparms["tb"]
170	318x	DTx <- function(x) {
171	636x	DTx.a <- (log(100/(100 - x)))/k1
172	636x	DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
173	339x	if (DTx.a < tb) DTx <- DTx.a
174	297x	else DTx <- DTx.b
175	636x	return(DTx)
176		}
177	318x	DT50 <- DTx(50)
178	318x	DT90 <- DTx(90)
179	318x	DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
180	318x	DT50_k1 = log(2)/k1
181	318x	DT50_k2 = log(2)/k2
182	318x	ep$distimes[obs_var, c("DT50back")] = DT50_back
183	318x	ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
184	318x	ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
185		}
186	73858x	if (type == "SFORB") {
187		# FOCUS kinetics (2006), p. 60 f
188	2616x	k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(degparms), value=TRUE)
189	2616x	k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
190	2616x	k_1output = sum(degparms[k_out_names])
191	2616x	k_12 = degparms[paste("k", obs_var, "free", "bound", sep="_")]
192	2616x	k_21 = degparms[paste("k", obs_var, "bound", "free", sep="_")]
193
194	2616x	sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21)
195	2616x	b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
196	2616x	b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
197	2616x	g = (k_12 + k_21 - b1)/(b2 - b1)
198
199	2616x	DT50_b1 = log(2)/b1
200	2616x	DT50_b2 = log(2)/b2
201	2616x	DT90_b1 = log(10)/b1
202	2616x	DT90_b2 = log(10)/b2
203
204	2616x	SFORB_fraction = function(t) {
205	60096x	g * exp(-b1 * t) + (1 - g) * exp(-b2 * t)
206		}
207
208	2616x	f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2
209	2616x	log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum,
210	2616x	silent = TRUE)
211	2616x	f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2
212	2616x	log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum,
213	2616x	silent = TRUE)
214
215	2616x	DT50 = if (inherits(log_DT50, "try-error")) NA
216	2616x	else exp(log_DT50)
217	2616x	DT90 = if (inherits(log_DT90, "try-error")) NA
218	2616x	else exp(log_DT90)
219
220	2616x	DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
221
222	2616x	for (k_out_name in k_out_names)
223		{
224	2618x	ep$ff[[sub("k_", "", k_out_name)]] = degparms[[k_out_name]] / k_1output
225		}
226
227		# Return the eigenvalues for comparison with DFOP rate constants
228	2616x	ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
229	2616x	ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
230		# Return g for comparison with DFOP
231	2616x	ep$SFORB[[paste(obs_var, "g", sep="_")]] = g
232
233	2616x	ep$distimes[obs_var, c("DT50back")] = DT50_back
234	2616x	ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1
235	2616x	ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2
236		}
237	73858x	if (type == "logistic") {
238		# FOCUS kinetics (2014) p. 67
239	153x	kmax = degparms["kmax"]
240	153x	k0 = degparms["k0"]
241	153x	r = degparms["r"]
242	153x	DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))
243	153x	DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))
244
245	153x	DT50_k0 = log(2)/k0
246	153x	DT50_kmax = log(2)/kmax
247	153x	ep$distimes[obs_var, c("DT50_k0")] = DT50_k0
248	153x	ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax
249		}
250	73858x	ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)
251		}
252	38846x	if (length(ep$ff) == 0) ep$ff <- NULL
253	53592x	if (length(ep$SFORB) == 0) ep$SFORB <- NULL
254	56208x	return(ep)
255		}

1		#' Fit nonlinear mixed-effects models built from one or more kinetic
2		#' degradation models and one or more error models
3		#'
4		#' The name of the methods expresses that (multiple) hierarchichal
5		#' (also known as multilevel) multicompartment kinetic models are
6		#' fitted. Our kinetic models are nonlinear, so we can use various nonlinear
7		#' mixed-effects model fitting functions.
8		#'
9		#' @param objects A list of [mmkin] objects containing fits of the same
10		#' degradation models to the same data, but using different error models.
11		#' Alternatively, a single [mmkin] object containing fits of several
12		#' degradation models to the same data
13		#' @param backend The backend to be used for fitting. Currently, only saemix is
14		#' supported
15		#' @param no_random_effect Default is NULL and will be passed to [saem]. If a
16		#' character vector is supplied, it will be passed to all calls to [saem],
17		#' which will exclude random effects for all matching parameters. Alternatively,
18		#' a list of character vectors or an object of class [illparms.mhmkin] can be
19		#' specified. They have to have the same dimensions that the return object of
20		#' the current call will have, i.e. the number of rows must match the number
21		#' of degradation models in the mmkin object(s), and the number of columns must
22		#' match the number of error models used in the mmkin object(s).
23		#' @param algorithm The algorithm to be used for fitting (currently not used)
24		#' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
25		#' model fitting function.
26		#' @param cores The number of cores to be used for multicore processing. This
27		#' is only used when the \code{cluster} argument is \code{NULL}. On Windows
28		#' machines, cores > 1 is not supported, you need to use the \code{cluster}
29		#' argument to use multiple logical processors. Per default, all cores detected
30		#' by [parallel::detectCores()] are used, except on Windows where the default
31		#' is 1.
32		#' @param cluster A cluster as returned by [makeCluster] to be used for
33		#' parallel execution.
34		#' @importFrom parallel mclapply parLapply detectCores
35		#' @return A two-dimensional [array] of fit objects and/or try-errors that can
36		#' be indexed using the degradation model names for the first index (row index)
37		#' and the error model names for the second index (column index), with class
38		#' attribute 'mhmkin'.
39		#' @author Johannes Ranke
40		#' @seealso \code{\link{[.mhmkin}} for subsetting [mhmkin] objects
41		#' @export
42		mhmkin <- function(objects, ...) {
43	375x	UseMethod("mhmkin")
44		}
45
46		#' @export
47		#' @rdname mhmkin
48		mhmkin.mmkin <- function(objects, ...) {
49	!	mhmkin(list(objects), ...)
50		}
51
52		#' @export
53		#' @rdname mhmkin
54		#' @examples
55		#' \dontrun{
56		#' # We start with separate evaluations of all the first six datasets with two
57		#' # degradation models and two error models
58		#' f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
59		#' f_sep_tc <- update(f_sep_const, error_model = "tc")
60		#' # The mhmkin function sets up hierarchical degradation models aka
61		#' # nonlinear mixed-effects models for all four combinations, specifying
62		#' # uncorrelated random effects for all degradation parameters
63		#' f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
64		#' status(f_saem_1)
65		#' # The 'illparms' function shows that in all hierarchical fits, at least
66		#' # one random effect is ill-defined (the confidence interval for the
67		#' # random effect expressed as standard deviation includes zero)
68		#' illparms(f_saem_1)
69		#' # Therefore we repeat the fits, excluding the ill-defined random effects
70		#' f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
71		#' status(f_saem_2)
72		#' illparms(f_saem_2)
73		#' # Model comparisons show that FOMC with two-component error is preferable,
74		#' # and confirms our reduction of the default parameter model
75		#' anova(f_saem_1)
76		#' anova(f_saem_2)
77		#' # The convergence plot for the selected model looks fine
78		#' saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
79		#' # The plot of predictions versus data shows that we have a pretty data-rich
80		#' # situation with homogeneous distribution of residuals, because we used the
81		#' # same degradation model, error model and parameter distribution model that
82		#' # was used in the data generation.
83		#' plot(f_saem_2[["FOMC", "tc"]])
84		#' # We can specify the same parameter model reductions manually
85		#' no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
86		#' dim(no_ranef) <- c(2, 2)
87		#' f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
88		#' anova(f_saem_2m)
89		#' }
90		mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
91		no_random_effect = NULL,
92		...,
93		cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
94		{
95	375x	call <- match.call()
96	375x	dot_args <- list(...)
97	375x	backend_function <- switch(backend,
98	375x	saemix = "saem"
99		)
100
101	375x	deg_models <- lapply(objects[[1]][, 1], function(x) x$mkinmod)
102	375x	names(deg_models) <- dimnames(objects[[1]])$model
103	375x	n.deg <- length(deg_models)
104
105	375x	ds <- lapply(objects[[1]][1, ], function(x) x$data)
106
107	375x	for (other in objects[-1]) {
108		# Check if the degradation models in all objects are the same
109	375x	for (deg_model_name in names(deg_models)) {
110	750x	if (!all.equal(other[[deg_model_name, 1]]$mkinmod$spec,
111	750x	deg_models[[deg_model_name]]$spec))
112		{
113	!	stop("The mmkin objects have to be based on the same degradation models")
114		}
115		}
116		# Check if they have been fitted to the same dataset
117	375x	other_object_ds <- lapply(other[1, ], function(x) x$data)
118	375x	for (i in 1:length(ds)) {
119	2250x	if (!all.equal(ds[[i]][c("time", "variable", "observed")],
120	2250x	other_object_ds[[i]][c("time", "variable", "observed")]))
121		{
122	!	stop("The mmkin objects have to be fitted to the same datasets")
123		}
124		}
125		}
126
127	375x	n.o <- length(objects)
128
129	375x	error_models = sapply(objects, function(x) x[[1]]$err_mod)
130	375x	n.e <- length(error_models)
131
132	375x	n.fits <- n.deg * n.e
133	375x	fit_indices <- matrix(1:n.fits, ncol = n.e)
134	375x	dimnames(fit_indices) <- list(degradation = names(deg_models),
135	375x	error = error_models)
136
137	375x	if (is.null(no_random_effect) \|\| is.null(dim(no_random_effect))) {
138	129x	no_ranef <- rep(list(no_random_effect), n.fits)
139	129x	dim(no_ranef) <- dim(fit_indices)
140		} else {
141	246x	if (!identical(dim(no_random_effect), dim(fit_indices))) {
142	!	stop("Dimensions of argument 'no_random_effect' are not suitable")
143		}
144	246x	if (is(no_random_effect, "illparms.mhmkin")) {
145	125x	no_ranef_dim <- dim(no_random_effect)
146	125x	no_ranef <- lapply(no_random_effect, function(x) {
147	500x	no_ranef_split <- strsplit(x, ", ")
148	500x	ret <- sapply(no_ranef_split, function(y) {
149	500x	gsub("sd\\((.*)\\)", "\\1", y)
150		})
151	500x	return(ret)
152		})
153	125x	dim(no_ranef) <- no_ranef_dim
154		} else {
155	121x	no_ranef <- no_random_effect
156		}
157		}
158
159	375x	fit_function <- function(fit_index) {
160	12x	w <- which(fit_indices == fit_index, arr.ind = TRUE)
161	12x	deg_index <- w[1]
162	12x	error_index <- w[2]
163	12x	mmkin_row <- objects[[error_index]][deg_index, ]
164	12x	res <- try(do.call(backend_function,
165	12x	args = c(
166	12x	list(mmkin_row),
167	12x	dot_args,
168	12x	list(no_random_effect = no_ranef[[deg_index, error_index]]))))
169	12x	return(res)
170		}
171
172
173	375x	fit_time <- system.time({
174	375x	if (is.null(cluster)) {
175	375x	results <- parallel::mclapply(as.list(1:n.fits), fit_function,
176	375x	mc.cores = cores, mc.preschedule = FALSE)
177		} else {
178	!	results <- parallel::parLapply(cluster, as.list(1:n.fits), fit_function)
179		}
180		})
181
182	363x	attributes(results) <- attributes(fit_indices)
183	363x	attr(results, "call") <- call
184	363x	attr(results, "time") <- fit_time
185	363x	class(results) <- switch(backend,
186	363x	saemix = c("mhmkin.saem.mmkin", "mhmkin")
187		)
188	363x	return(results)
189		}
190
191		#' Subsetting method for mhmkin objects
192		#'
193		#' @param x An [mhmkin] object.
194		#' @param i Row index selecting the fits for specific models
195		#' @param j Column index selecting the fits to specific datasets
196		#' @param drop If FALSE, the method always returns an mhmkin object, otherwise
197		#' either a list of fit objects or a single fit object.
198		#' @return An object inheriting from \code{\link{mhmkin}}.
199		#' @rdname mhmkin
200		#' @export
201		`[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
202	!	original_class <- class(x)
203	!	class(x) <- NULL
204	!	x_sub <- x[i, j, drop = drop]
205
206	!	if (!drop) {
207	!	class(x_sub) <- original_class
208		}
209	!	return(x_sub)
210		}
211
212		#' Print method for mhmkin objects
213		#'
214		#' @rdname mhmkin
215		#' @export
216		print.mhmkin <- function(x, ...) {
217	125x	cat("<mhmkin> object\n")
218	125x	cat("Status of individual fits:\n\n")
219	125x	print(status(x))
220		}
221
222		#' Check if fit within an mhmkin object failed
223		#' @param x The object to be checked
224		check_failed <- function(x) {
225	1936x	if (inherits(x, "try-error")) {
226	!	return(TRUE)
227		} else {
228	1936x	if (inherits(x$so, "try-error")) {
229	!	return(TRUE)
230		} else {
231	1936x	return(FALSE)
232		}
233		}
234		}
235
236		#' @export
237		AIC.mhmkin <- function(object, ..., k = 2) {
238	125x	res <- sapply(object, function(x) {
239	!	if (check_failed(x)) return(NA)
240	500x	else return(AIC(x$so, k = k))
241		})
242	125x	dim(res) <- dim(object)
243	125x	dimnames(res) <- dimnames(object)
244	125x	return(res)
245		}
246
247		#' @export
248		BIC.mhmkin <- function(object, ...) {
249	125x	res <- sapply(object, function(x) {
250	!	if (check_failed(x)) return(NA)
251	500x	else return(BIC(x$so))
252		})
253	125x	dim(res) <- dim(object)
254	125x	dimnames(res) <- dimnames(object)
255	125x	return(res)
256		}
257
258		#' @export
259		update.mhmkin <- function(object, ..., evaluate = TRUE) {
260	246x	call <- attr(object, "call")
261		# For some reason we get mhkin.list in call[[1]] when using mhmkin from the
262		# loaded package so we need to fix this so we do not have to export
263		# mhmkin.list in addition to the S3 method mhmkin
264	246x	call[[1]] <- mhmkin
265
266	246x	update_arguments <- match.call(expand.dots = FALSE)$...
267
268	246x	if (length(update_arguments) > 0) {
269	246x	update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
270		}
271
272	246x	for (a in names(update_arguments)[update_arguments_in_call]) {
273	!	call[[a]] <- update_arguments[[a]]
274		}
275
276	246x	update_arguments_not_in_call <- !update_arguments_in_call
277	246x	if(any(update_arguments_not_in_call)) {
278	246x	call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
279	246x	call <- as.call(call)
280		}
281	246x	if(evaluate) eval(call, parent.frame())
282	!	else call
283		}
284
285		#' @export
286		anova.mhmkin <- function(object, ...,
287		method = c("is", "lin", "gq"), test = FALSE, model.names = "auto") {
288	234x	if (identical(model.names, "auto")) {
289	234x	model.names <- outer(rownames(object), colnames(object), paste)
290		}
291	234x	failed_index <- which(sapply(object, check_failed), arr.ind = TRUE)
292	234x	if (length(failed_index > 0)) {
293	!	rlang::inject(anova(!!!(object[-failed_index]), method = method, test = test,
294	!	model.names = model.names[-failed_index]))
295		} else {
296	234x	rlang::inject(anova(!!!(object), method = method, test = test,
297	234x	model.names = model.names))
298		}
299		}
300

1		#' Evaluate parent kinetics using the NAFTA guidance
2		#'
3		#' The function fits the SFO, IORE and DFOP models using \code{\link{mmkin}}
4		#' and returns an object of class \code{nafta} that has methods for printing
5		#' and plotting.
6		#'
7		#' @param ds A dataframe that must contain one variable called "time" with the
8		#' time values specified by the \code{time} argument, one column called
9		#' "name" with the grouping of the observed values, and finally one column of
10		#' observed values called "value".
11		#' @param title Optional title of the dataset
12		#' @param quiet Should the evaluation text be shown?
13		#' @param \dots Further arguments passed to \code{\link{mmkin}} (not for the
14		#' printing method).
15		#' @importFrom stats qf
16		#' @return An list of class \code{nafta}. The list element named "mmkin" is the
17		#' \code{\link{mmkin}} object containing the fits of the three models. The
18		#' list element named "title" contains the title of the dataset used. The
19		#' list element "data" contains the dataset used in the fits.
20		#' @author Johannes Ranke
21		#' @source NAFTA (2011) Guidance for evaluating and calculating degradation
22		#' kinetics in environmental media. NAFTA Technical Working Group on
23		#' Pesticides
24		#' \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/guidance-evaluating-and-calculating-degradation}
25		#' accessed 2019-02-22
26		#'
27		#' US EPA (2015) Standard Operating Procedure for Using the NAFTA Guidance to
28		#' Calculate Representative Half-life Values and Characterizing Pesticide
29		#' Degradation
30		#' \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/standard-operating-procedure-using-nafta-guidance}
31		#' @examples
32		#'
33		#' nafta_evaluation <- nafta(NAFTA_SOP_Appendix_D, cores = 1)
34		#' print(nafta_evaluation)
35		#' plot(nafta_evaluation)
36		#'
37		#' @export
38		nafta <- function(ds, title = NA, quiet = FALSE, ...) {
39	264x	if (length(levels(ds$name)) > 1) {
40	88x	stop("The NAFTA procedure is only defined for decline data for a single compound")
41		}
42	176x	n <- nrow(subset(ds, !is.na(value)))
43	176x	models <- c("SFO", "IORE", "DFOP")
44
45	176x	result <- list(title = title, data = ds)
46	176x	result$mmkin <- mmkin(models, list(ds), quiet = TRUE, ...)
47
48	176x	distimes <- lapply(result$mmkin, function(x) as.numeric(endpoints(x)$distimes["parent", ]))
49
50	176x	result$distimes <- matrix(NA, nrow = 3, ncol = 3,
51	176x	dimnames = list(models, c("DT50", "DT90", "DT50_rep")))
52	176x	result$distimes["SFO", ] <- distimes[[1]][c(1, 2, 1)]
53	176x	result$distimes["IORE", ] <- distimes[[2]][c(1, 2, 3)]
54	176x	result$distimes["DFOP", ] <- distimes[[3]][c(1, 2, 5)]
55
56		# Get parameters with statistics
57	176x	result$parameters <- lapply(result$mmkin, function(x) {
58	528x	summary(x)$bpar[, c(1, 4:6)]
59		})
60	176x	names(result$parameters) <- models
61
62		# Compare the sum of squared residuals (SSR) to the upper bound of the
63		# confidence region of the SSR for the IORE model
64	176x	result$S <- sapply(result$mmkin, function(x) sum(x$data$residual^2))
65	176x	names(result$S) <- c("SFO", "IORE", "DFOP")
66		# Equation (3) on p. 3
67	176x	p <- 3
68	176x	result$S["IORE"]
69	176x	result$S_c <- result$S[["IORE"]] * (1 + p/(n - p) * qf(0.5, p, n - p))
70
71	176x	result$t_rep <- .evaluate_nafta_results(result$S, result$S_c,
72	176x	result$distimes, quiet = quiet)
73
74	176x	class(result) <- "nafta"
75	176x	return(result)
76		}
77
78		#' Plot the results of the three models used in the NAFTA scheme.
79		#'
80		#' The plots are ordered with increasing complexity of the model in this
81		#' function (SFO, then IORE, then DFOP).
82		#'
83		#' Calls \code{\link{plot.mmkin}}.
84		#'
85		#' @param x An object of class \code{\link{nafta}}.
86		#' @param legend Should a legend be added?
87		#' @param main Possibility to override the main title of the plot.
88		#' @param \dots Further arguments passed to \code{\link{plot.mmkin}}.
89		#' @return The function is called for its side effect.
90		#' @author Johannes Ranke
91		#' @export
92		plot.nafta <- function(x, legend = FALSE, main = "auto", ...) {
93	176x	if (main == "auto") {
94	!	if (is.na(x$title)) main = ""
95	176x	else main = x$title
96		}
97	176x	plot(x$mmkin, ..., legend = legend, main = main)
98		}
99
100		#' Print nafta objects
101		#'
102		#' Print nafta objects. The results for the three models are printed in the
103		#' order of increasing model complexity, i.e. SFO, then IORE, and finally DFOP.
104		#'
105		#' @param x An \code{\link{nafta}} object.
106		#' @param digits Number of digits to be used for printing parameters and
107		#' dissipation times.
108		#' @rdname nafta
109		#' @export
110		print.nafta <- function(x, quiet = TRUE, digits = 3, ...) {
111	176x	cat("Sums of squares:\n")
112	176x	print(x$S)
113	176x	cat("\nCritical sum of squares for checking the SFO model:\n")
114	176x	print(x$S_c)
115	176x	cat("\nParameters:\n")
116	176x	print(x$parameters, digits = digits)
117	176x	t_rep <- .evaluate_nafta_results(x$S, x$S_c, x$distimes, quiet = quiet)
118	176x	cat("\nDTx values:\n")
119	176x	print(signif(x$distimes, digits = digits))
120	176x	cat("\nRepresentative half-life:\n")
121	176x	print(round(t_rep, 2))
122		}
123
124		.evaluate_nafta_results <- function(S, S_c, distimes, quiet = FALSE) {
125	352x	t_SFO <- distimes["IORE", "DT50"]
126	352x	t_IORE <- distimes["IORE", "DT50_rep"]
127	352x	t_DFOP2 <- distimes["DFOP", "DT50_rep"]
128
129	352x	if (S["SFO"] < S_c) {
130	!	if (!quiet) {
131	!	message("S_SFO is lower than the critical value S_c, use the SFO model")
132		}
133	!	t_rep <- t_SFO
134		} else {
135	352x	if (!quiet) {
136	88x	message("The SFO model is rejected as S_SFO is equal or higher than the critical value S_c")
137		}
138	352x	if (t_IORE < t_DFOP2) {
139	176x	if (!quiet) {
140	88x	message("The half-life obtained from the IORE model may be used")
141		}
142	176x	t_rep <- t_IORE
143		} else {
144	176x	if (!quiet) {
145	!	message("The representative half-life of the IORE model is longer than the one corresponding")
146	!	message("to the terminal degradation rate found with the DFOP model.")
147	!	message("The representative half-life obtained from the DFOP model may be used")
148		}
149	176x	t_rep <- t_DFOP2
150		}
151		}
152	352x	return(t_rep)
153		}

1		#' Summary method for class "nlme.mmkin"
2		#'
3		#' Lists model equations, initial parameter values, optimised parameters
4		#' for fixed effects (population), random effects (deviations from the
5		#' population mean) and residual error model, as well as the resulting
6		#' endpoints such as formation fractions and DT50 values. Optionally
7		#' (default is FALSE), the data are listed in full.
8		#'
9		#' @param object an object of class [nlme.mmkin]
10		#' @param x an object of class [summary.nlme.mmkin]
11		#' @param data logical, indicating whether the full data should be included in
12		#' the summary.
13		#' @param verbose Should the summary be verbose?
14		#' @param distimes logical, indicating whether DT50 and DT90 values should be
15		#' included.
16		#' @param alpha error level for confidence interval estimation from the t
17		#' distribution
18		#' @param digits Number of digits to use for printing
19		#' @param \dots optional arguments passed to methods like \code{print}.
20		#' @return The summary function returns a list based on the [nlme] object
21		#' obtained in the fit, with at least the following additional components
22		#' \item{nlmeversion, mkinversion, Rversion}{The nlme, mkin and R versions used}
23		#' \item{date.fit, date.summary}{The dates where the fit and the summary were
24		#' produced}
25		#' \item{diffs}{The differential equations used in the degradation model}
26		#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
27		#' \item{data}{The data}
28		#' \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals}
29		#' \item{confint_back}{Backtransformed parameters, with confidence intervals if available}
30		#' \item{ff}{The estimated formation fractions derived from the fitted
31		#' model.}
32		#' \item{distimes}{The DT50 and DT90 values for each observed variable.}
33		#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
34		#' The print method is called for its side effect, i.e. printing the summary.
35		#' @importFrom stats predict
36		#' @author Johannes Ranke for the mkin specific parts
37		#' José Pinheiro and Douglas Bates for the components inherited from nlme
38		#' @examples
39		#'
40		#' # Generate five datasets following SFO kinetics
41		#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
42		#' dt50_sfo_in_pop <- 50
43		#' k_in_pop <- log(2) / dt50_sfo_in_pop
44		#' set.seed(1234)
45		#' k_in <- rlnorm(5, log(k_in_pop), 0.5)
46		#' SFO <- mkinmod(parent = mkinsub("SFO"))
47		#'
48		#' pred_sfo <- function(k) {
49		#' mkinpredict(SFO,
50		#' c(k_parent = k),
51		#' c(parent = 100),
52		#' sampling_times)
53		#' }
54		#'
55		#' ds_sfo_mean <- lapply(k_in, pred_sfo)
56		#' names(ds_sfo_mean) <- paste("ds", 1:5)
57		#'
58		#' set.seed(12345)
59		#' ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) {
60		#' add_err(ds,
61		#' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
62		#' n = 1)[[1]]
63		#' })
64		#'
65		#' \dontrun{
66		#' # Evaluate using mmkin and nlme
67		#' library(nlme)
68		#' f_mmkin <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1)
69		#' f_nlme <- nlme(f_mmkin)
70		#' summary(f_nlme, data = TRUE)
71		#' }
72		#'
73		#' @export
74		summary.nlme.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, alpha = 0.05, ...) {
75
76	319x	mod_vars <- names(object$mkinmod$diffs)
77
78	319x	confint_trans <- intervals(object, which = "fixed", level = 1 - alpha)$fixed
79	319x	attr(confint_trans, "label") <- NULL
80	319x	pnames <- rownames(confint_trans)
81
82	319x	bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
83	319x	object$transform_rates, object$transform_fractions)
84	319x	bpnames <- names(bp)
85
86		# variance-covariance estimates for fixed effects (from summary.lme)
87	319x	fixed <- fixef(object)
88	319x	stdFixed <- sqrt(diag(as.matrix(object$varFix)))
89	319x	object$corFixed <- array(
90	319x	t(object$varFix/stdFixed)/stdFixed,
91	319x	dim(object$varFix),
92	319x	list(names(fixed), names(fixed)))
93
94		# Transform boundaries of CI for one parameter at a time,
95		# with the exception of sets of formation fractions (single fractions are OK).
96	319x	f_names_skip <- character(0)
97	319x	for (box in mod_vars) { # Figure out sets of fractions to skip
98	436x	f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
99	436x	n_paths <- length(f_names)
100	!	if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
101		}
102
103	319x	confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
104	319x	dimnames = list(bpnames, colnames(confint_trans)))
105	319x	confint_back[, "est."] <- bp
106
107	319x	for (pname in pnames) {
108	1410x	if (!pname %in% f_names_skip) {
109	1410x	par.lower <- confint_trans[pname, "lower"]
110	1410x	par.upper <- confint_trans[pname, "upper"]
111	1410x	names(par.lower) <- names(par.upper) <- pname
112	1410x	bpl <- backtransform_odeparms(par.lower, object$mkinmod,
113	1410x	object$transform_rates,
114	1410x	object$transform_fractions)
115	1410x	bpu <- backtransform_odeparms(par.upper, object$mkinmod,
116	1410x	object$transform_rates,
117	1410x	object$transform_fractions)
118	1410x	confint_back[names(bpl), "lower"] <- bpl
119	1410x	confint_back[names(bpu), "upper"] <- bpu
120		}
121		}
122
123	319x	object$confint_trans <- confint_trans
124	319x	object$confint_back <- confint_back
125
126	319x	object$date.summary = date()
127	319x	object$use_of_ff = object$mkinmod$use_of_ff
128	319x	object$error_model_algorithm = object$mmkin[[1]]$error_model_algorithm
129	319x	err_mod = object$mmkin[[1]]$err_mod
130
131	319x	object$diffs <- object$mkinmod$diffs
132	319x	object$print_data <- data
133	319x	object$data[["observed"]] <- object$data[["value"]]
134	319x	object$data[["value"]] <- NULL
135	319x	object$data[["predicted"]] <- predict(object)
136	319x	object$data[["residual"]] <- residuals(object, type = "response")
137	319x	if (is.null(object$modelStruct$varStruct)) {
138	!	object$data[["std"]] <- object$sigma
139		} else {
140	319x	object$data[["std"]] <- 1/attr(object$modelStruct$varStruct, "weights")
141		}
142	319x	object$data[["standardized"]] <- residuals(object, type = "pearson")
143	319x	object$verbose <- verbose
144
145	319x	object$fixed <- object$mmkin[[1]]$fixed
146	319x	object$AIC = AIC(object)
147	319x	object$BIC = BIC(object)
148	319x	object$logLik = logLik(object)
149
150	319x	ep <- endpoints(object)
151	319x	if (length(ep$ff) != 0)
152	117x	object$ff <- ep$ff
153	319x	if (distimes) object$distimes <- ep$distimes
154	!	if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB
155	319x	class(object) <- c("summary.nlme.mmkin", "nlme.mmkin", "nlme", "lme")
156	319x	return(object)
157		}
158
159		#' @rdname summary.nlme.mmkin
160		#' @export
161		print.summary.nlme.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
162	117x	cat("nlme version used for fitting: ", x$nlmeversion, "\n")
163	117x	cat("mkin version used for pre-fitting: ", x$mkinversion, "\n")
164	117x	cat("R version used for fitting: ", x$Rversion, "\n")
165
166	117x	cat("Date of fit: ", x$date.fit, "\n")
167	117x	cat("Date of summary:", x$date.summary, "\n")
168
169	117x	cat("\nEquations:\n")
170	117x	nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
171	117x	writeLines(strwrap(nice_diffs, exdent = 11))
172
173	117x	cat("\nData:\n")
174	117x	cat(nrow(x$data), "observations of",
175	117x	length(unique(x$data$name)), "variable(s) grouped in",
176	117x	length(unique(x$data$ds)), "datasets\n")
177
178	117x	cat("\nModel predictions using solution type", x$solution_type, "\n")
179
180	117x	cat("\nFitted in", x$time[["elapsed"]], "s using", x$numIter, "iterations\n")
181
182	117x	cat("\nVariance model: ")
183	117x	cat(switch(x$err_mod,
184	117x	const = "Constant variance",
185	117x	obs = "Variance unique to each observed variable",
186	117x	tc = "Two-component variance function"), "\n")
187
188	117x	cat("\nMean of starting values for individual parameters:\n")
189	117x	print(x$mean_dp_start, digits = digits)
190
191	117x	cat("\nFixed degradation parameter values:\n")
192	117x	if(length(x$fixed$value) == 0) cat("None\n")
193	!	else print(x$fixed, digits = digits)
194
195	117x	cat("\nResults:\n\n")
196	117x	print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
197	117x	row.names = " "), digits = digits, ...)
198
199	117x	cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
200	117x	print(x$confint_trans, digits = digits, ...)
201
202	117x	if (nrow(x$confint_trans) > 1) {
203	117x	corr <- x$corFixed
204	117x	class(corr) <- "correlation"
205	117x	print(corr, title = "\nCorrelation:", rdig = digits, ...)
206		}
207
208	117x	cat("\n") # Random effects
209	117x	print(summary(x$modelStruct), sigma = x$sigma,
210	117x	reEstimates = x$coef$random, digits = digits, verbose = verbose, ...)
211
212	117x	cat("\nBacktransformed parameters with asymmetric confidence intervals:\n")
213	117x	print(x$confint_back, digits = digits, ...)
214
215
216	117x	printSFORB <- !is.null(x$SFORB)
217	117x	if(printSFORB){
218	!	cat("\nEstimated Eigenvalues of SFORB model(s):\n")
219	!	print(x$SFORB, digits = digits,...)
220		}
221
222	117x	printff <- !is.null(x$ff)
223	117x	if(printff){
224	!	cat("\nResulting formation fractions:\n")
225	!	print(data.frame(ff = x$ff), digits = digits, ...)
226		}
227
228	117x	printdistimes <- !is.null(x$distimes)
229	117x	if(printdistimes){
230	117x	cat("\nEstimated disappearance times:\n")
231	117x	print(x$distimes, digits = digits, ...)
232		}
233
234	117x	if (x$print_data){
235	!	cat("\nData:\n")
236	!	print(format(x$data, digits = digits, ...), row.names = FALSE)
237		}
238
239	117x	invisible(x)
240		}

1		utils::globalVariables(c("variable", "residual"))
2
3		#' Function to plot squared residuals and the error model for an mkin object
4		#'
5		#' This function plots the squared residuals for the specified subset of the
6		#' observed variables from an mkinfit object. In addition, one or more dashed
7		#' line(s) show the fitted error model. A combined plot of the fitted model
8		#' and this error model plot can be obtained with \code{\link{plot.mkinfit}}
9		#' using the argument \code{show_errplot = TRUE}.
10		#'
11		#' @param object A fit represented in an \code{\link{mkinfit}} object.
12		#' @param obs_vars A character vector of names of the observed variables for
13		#' which residuals should be plotted. Defaults to all observed variables in
14		#' the model
15		#' @param xlim plot range in x direction.
16		#' @param xlab Label for the x axis.
17		#' @param ylab Label for the y axis.
18		#' @param maxy Maximum value of the residuals. This is used for the scaling of
19		#' the y axis and defaults to "auto".
20		#' @param legend Should a legend be plotted?
21		#' @param lpos Where should the legend be placed? Default is "topright". Will
22		#' be passed on to \code{\link{legend}}.
23		#' @param col_obs Colors for the observed variables.
24		#' @param pch_obs Symbols to be used for the observed variables.
25		#' @param frame Should a frame be drawn around the plots?
26		#' @param \dots further arguments passed to \code{\link{plot}}.
27		#' @return Nothing is returned by this function, as it is called for its side
28		#' effect, namely to produce a plot.
29		#' @author Johannes Ranke
30		#' @seealso \code{\link{mkinplot}}, for a way to plot the data and the fitted
31		#' lines of the mkinfit object.
32		#' @keywords hplot
33		#' @examples
34		#'
35		#' \dontrun{
36		#' model <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
37		#' fit <- mkinfit(model, FOCUS_2006_D, error_model = "tc", quiet = TRUE)
38		#' mkinerrplot(fit)
39		#' }
40		#'
41		#' @export
42		mkinerrplot <- function (object,
43		obs_vars = names(object$mkinmod$map),
44		xlim = c(0, 1.1 * max(object$data$predicted)),
45		xlab = "Predicted", ylab = "Squared residual",
46		maxy = "auto", legend= TRUE, lpos = "topright",
47		col_obs = "auto", pch_obs = "auto",
48		frame = TRUE,
49		...)
50		{
51	275x	obs_vars_all <- as.character(unique(object$data$variable))
52
53	275x	if (length(obs_vars) > 0){
54	275x	obs_vars <- intersect(obs_vars_all, obs_vars)
55	!	} else obs_vars <- obs_vars_all
56
57	275x	residuals <- subset(object$data, variable %in% obs_vars, residual)
58
59	275x	if (maxy == "auto") maxy = max(residuals^2, na.rm = TRUE)
60
61		# Set colors and symbols
62	275x	if (col_obs[1] == "auto") {
63	70x	col_obs <- 1:length(obs_vars)
64		}
65
66	275x	if (pch_obs[1] == "auto") {
67	70x	pch_obs <- 1:length(obs_vars)
68		}
69	275x	names(col_obs) <- names(pch_obs) <- obs_vars
70
71	275x	plot(0, type = "n",
72	275x	xlab = xlab, ylab = ylab,
73	275x	xlim = xlim,
74	275x	ylim = c(0, 1.2 * maxy), frame = frame, ...)
75
76	275x	for(obs_var in obs_vars){
77	410x	residuals_plot <- subset(object$data, variable == obs_var, c("predicted", "residual"))
78	410x	points(residuals_plot[["predicted"]],
79	410x	residuals_plot[["residual"]]^2,
80	410x	pch = pch_obs[obs_var], col = col_obs[obs_var])
81		}
82
83	275x	if (object$err_mod == "const") {
84	140x	abline(h = object$errparms^2, lty = 2, col = 1)
85		}
86	275x	if (object$err_mod == "obs") {
87	65x	for (obs_var in obs_vars) {
88	130x	sigma_name = paste0("sigma_", obs_var)
89	130x	abline(h = object$errparms[sigma_name]^2, lty = 2,
90	130x	col = col_obs[obs_var])
91		}
92		}
93	275x	if (object$err_mod == "tc") {
94	70x	sigma_plot <- function(predicted) {
95	70x	sigma_twocomp(predicted,
96	70x	sigma_low = object$errparms[1],
97	70x	rsd_high = object$errparms[2])^2
98		}
99	70x	plot(sigma_plot, from = 0, to = max(object$data$predicted),
100	70x	add = TRUE, lty = 2, col = 1)
101		}
102
103	275x	if (legend == TRUE) {
104	70x	legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
105	70x	col = col_obs[obs_vars], pch = pch_obs[obs_vars])
106		}
107		}

1		utils::globalVariables("ds")
2
3		#' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
4		#'
5		#' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin]
6		#' @param i A numeric index to select datasets for which to plot the individual predictions,
7		#' in case plots get too large
8		#' @inheritParams plot.mkinfit
9		#' @param standardized Should the residuals be standardized? Only takes effect if
10		#' `resplot = "time"`.
11		#' @param pop_curves Per default, one population curve is drawn in case
12		#' population parameters are fitted by the model, e.g. for saem objects.
13		#' In case there is a covariate model, the behaviour depends on the value
14		#' of 'covariates'
15		#' @param covariates Data frame with covariate values for all variables in
16		#' any covariate models in the object. If given, it overrides 'covariate_quantiles'.
17		#' Each line in the data frame will result in a line drawn for the population.
18		#' Rownames are used in the legend to label the lines.
19		#' @param covariate_quantiles This argument only has an effect if the fitted
20		#' object has covariate models. If so, the default is to show three population
21		#' curves, for the 5th percentile, the 50th percentile and the 95th percentile
22		#' of the covariate values used for fitting the model.
23		#' @note Covariate models are currently only supported for saem.mmkin objects.
24		#' @param pred_over Named list of alternative predictions as obtained
25		#' from [mkinpredict] with a compatible [mkinmod].
26		#' @param test_log_parms Passed to [mean_degparms] in the case of an
27		#' [mixed.mmkin] object
28		#' @param conf.level Passed to [mean_degparms] in the case of an
29		#' [mixed.mmkin] object
30		#' @param default_log_parms Passed to [mean_degparms] in the case of an
31		#' [mixed.mmkin] object
32		#' @param rel.height.legend The relative height of the legend shown on top
33		#' @param rel.height.bottom The relative height of the bottom plot row
34		#' @param ymax Vector of maximum y axis values
35		#' @param ncol.legend Number of columns to use in the legend
36		#' @param nrow.legend Number of rows to use in the legend
37		#' @param resplot Should the residuals plotted against time or against
38		#' predicted values?
39		#' @param col_ds Colors used for plotting the observed data and the
40		#' corresponding model prediction lines for the different datasets.
41		#' @param pch_ds Symbols to be used for plotting the data.
42		#' @param lty_ds Line types to be used for the model predictions.
43		#' @importFrom stats coefficients
44		#' @return The function is called for its side effect.
45		#' @author Johannes Ranke
46		#' @examples
47		#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
48		#' function(x) x$data[c("name", "time", "value")])
49		#' names(ds) <- paste0("ds ", 6:10)
50		#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
51		#' A1 = mkinsub("SFO"), quiet = TRUE)
52		#' \dontrun{
53		#' f <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE)
54		#' plot(f[, 3:4], standardized = TRUE)
55		#'
56		#' # For this fit we need to increase pnlsMaxiter, and we increase the
57		#' # tolerance in order to speed up the fit for this example evaluation
58		#' # It still takes 20 seconds to run
59		#' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
60		#' plot(f_nlme)
61		#'
62		#' f_saem <- saem(f, transformations = "saemix")
63		#' plot(f_saem)
64		#'
65		#' f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs")
66		#' f_nlmix <- nlmix(f_obs)
67		#' plot(f_nlmix)
68		#'
69		#' # We can overlay the two variants if we generate predictions
70		#' pred_nlme <- mkinpredict(dfop_sfo,
71		#' f_nlme$bparms.optim[-1],
72		#' c(parent = f_nlme$bparms.optim[[1]], A1 = 0),
73		#' seq(0, 180, by = 0.2))
74		#' plot(f_saem, pred_over = list(nlme = pred_nlme))
75		#' }
76		#' @export
77		plot.mixed.mmkin <- function(x,
78		i = 1:ncol(x$mmkin),
79		obs_vars = names(x$mkinmod$map),
80		standardized = TRUE,
81		covariates = NULL,
82		covariate_quantiles = c(0.5, 0.05, 0.95),
83		xlab = "Time",
84		xlim = range(x$data$time),
85		resplot = c("predicted", "time"),
86		pop_curves = "auto",
87		pred_over = NULL,
88		test_log_parms = FALSE,
89		conf.level = 0.6,
90		default_log_parms = NA,
91		ymax = "auto", maxabs = "auto",
92		ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
93		nrow.legend = ceiling((length(i) + 1) / ncol.legend),
94		rel.height.legend = 0.02 + 0.07 * nrow.legend,
95		rel.height.bottom = 1.1,
96		pch_ds = 1:length(i),
97		col_ds = pch_ds + 1,
98		lty_ds = col_ds,
99		frame = TRUE, ...
100		)
101		{
102		# Prepare parameters and data
103	283x	fit_1 <- x$mmkin[[1]]
104	283x	ds_names <- colnames(x$mmkin)
105
106	283x	backtransform = TRUE
107
108	283x	if (identical(class(x), "mixed.mmkin")) {
109	65x	if (identical(pop_curves, "auto")) {
110	!	pop_curves <- FALSE
111		} else {
112	65x	pop_curves <- TRUE
113		}
114	65x	if (pop_curves) {
115	65x	degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms,
116	65x	conf.level = conf.level, default_log_parms = default_log_parms)
117		}
118
119	65x	degparms_tmp <- parms(x$mmkin, transformed = TRUE)
120	65x	degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ]))
121	65x	residual_type = ifelse(standardized, "standardized", "residual")
122	65x	residuals <- x$data[[residual_type]]
123		}
124
125	283x	if (inherits(x, "nlme.mmkin")) {
126	65x	if (identical(pop_curves, "auto")) {
127	65x	pop_curves <- TRUE
128		} else {
129	!	pop_curves <- FALSE
130		}
131	65x	degparms_i <- coefficients(x)
132	65x	degparms_pop <- nlme::fixef(x)
133	65x	residuals <- residuals(x,
134	65x	type = ifelse(standardized, "pearson", "response"))
135		}
136
137	283x	if (inherits(x, "saem.mmkin")) {
138	65x	if (x$transformations == "saemix") backtransform = FALSE
139	153x	psi <- saemix::psi(x$so)
140	153x	rownames(psi) <- x$saemix_ds_order
141	153x	degparms_i <- psi[ds_names, ]
142	153x	degparms_i_names <- colnames(degparms_i)
143	153x	residual_type = ifelse(standardized, "standardized", "residual")
144	153x	residuals <- x$data[[residual_type]]
145
146	153x	if (identical(pop_curves, "auto")) {
147	153x	if (length(x$covariate_models) == 0) {
148	153x	degparms_pop <- x$so@results@fixed.effects
149	153x	names(degparms_pop) <- degparms_i_names
150	153x	pop_curves <- TRUE
151		} else {
152	!	if (is.null(covariates)) {
153	!	covariates = as.data.frame(
154	!	apply(x$covariates, 2, quantile,
155	!	covariate_quantiles, simplify = FALSE))
156	!	rownames(covariates) <- paste(
157	!	ifelse(length(x$covariate_models) == 1,
158	!	"Covariate", "Covariates"),
159	!	rownames(covariates))
160		}
161	!	degparms_pop <- parms(x, covariates = covariates)
162	!	pop_curves <- TRUE
163		}
164		} else {
165	!	pop_curves <- FALSE
166		}
167		}
168
169	283x	if (pop_curves) {
170		# Make sure degparms_pop is a matrix, columns corresponding to population curve(s)
171	283x	if (is.null(dim(degparms_pop))) {
172	283x	degparms_pop <- matrix(degparms_pop, ncol = 1,
173	283x	dimnames = list(names(degparms_pop), "Population"))
174		}
175		}
176
177	283x	degparms_fixed <- fit_1$fixed$value
178	283x	names(degparms_fixed) <- rownames(fit_1$fixed)
179	283x	degparms_all <- cbind(as.matrix(degparms_i),
180	283x	matrix(rep(degparms_fixed, nrow(degparms_i)),
181	283x	ncol = length(degparms_fixed),
182	283x	nrow = nrow(degparms_i), byrow = TRUE))
183	283x	degparms_all_names <- c(names(degparms_i), names(degparms_fixed))
184	283x	colnames(degparms_all) <- degparms_all_names
185
186	283x	odeini_names <- grep("_0$", degparms_all_names, value = TRUE)
187	283x	odeparms_names <- setdiff(degparms_all_names, odeini_names)
188
189	283x	observed <- cbind(x$data[c("ds", "name", "time", "value")],
190	283x	residual = residuals)
191
192	283x	solution_type = fit_1$solution_type
193
194	283x	outtimes <- sort(unique(c(x$data$time,
195	283x	seq(xlim[1], xlim[2], length.out = 50))))
196
197	283x	pred_list <- lapply(i, function(ds_i) {
198	2945x	odeparms_trans <- degparms_all[ds_i, odeparms_names]
199	2945x	names(odeparms_trans) <- odeparms_names # needed if only one odeparm
200	2945x	if (backtransform) {
201	2620x	odeparms <- backtransform_odeparms(odeparms_trans,
202	2620x	x$mkinmod,
203	2620x	transform_rates = fit_1$transform_rates,
204	2620x	transform_fractions = fit_1$transform_fractions)
205		} else {
206	325x	odeparms <- odeparms_trans
207		}
208
209	2945x	odeini <- degparms_all[ds_i, odeini_names]
210	2945x	names(odeini) <- gsub("_0", "", odeini_names)
211
212	2945x	out <- mkinpredict(x$mkinmod, odeparms, odeini,
213	2945x	outtimes, solution_type = solution_type,
214	2945x	atol = fit_1$atol, rtol = fit_1$rtol)
215		})
216	283x	names(pred_list) <- ds_names[i]
217	283x	pred_ds <- vctrs::vec_rbind(!!!pred_list, .names_to = "ds")
218
219	283x	if (pop_curves) {
220	283x	pred_list_pop <- lapply(1:ncol(degparms_pop), function(cov_i) {
221	283x	degparms_all_pop_i <- c(degparms_pop[, cov_i], degparms_fixed)
222	283x	odeparms_pop_trans_i <- degparms_all_pop_i[odeparms_names]
223	283x	names(odeparms_pop_trans_i) <- odeparms_names # needed if only one odeparm
224	283x	if (backtransform) {
225	218x	odeparms_pop_i <- backtransform_odeparms(odeparms_pop_trans_i,
226	218x	x$mkinmod,
227	218x	transform_rates = fit_1$transform_rates,
228	218x	transform_fractions = fit_1$transform_fractions)
229		} else {
230	65x	odeparms_pop_i <- odeparms_pop_trans_i
231		}
232
233	283x	odeini <- degparms_all_pop_i[odeini_names]
234	283x	names(odeini) <- gsub("_0", "", odeini_names)
235
236	283x	out <- mkinpredict(x$mkinmod, odeparms_pop_i, odeini,
237	283x	outtimes, solution_type = solution_type,
238	283x	atol = fit_1$atol, rtol = fit_1$rtol)
239		})
240	283x	names(pred_list_pop) <- colnames(degparms_pop)
241
242		} else {
243	!	pred_list_pop <- NULL
244		}
245
246		# Start of graphical section
247	283x	oldpar <- par(no.readonly = TRUE)
248	283x	on.exit(par(oldpar, no.readonly = TRUE))
249
250	283x	n_plot_rows = length(obs_vars)
251	283x	n_plots = n_plot_rows * 2
252
253		# Set relative plot heights, so the first plot row is the norm
254	283x	rel.heights <- if (n_plot_rows > 1) {
255	218x	c(rel.height.legend, c(rep(1, n_plot_rows - 1), rel.height.bottom))
256		} else {
257	65x	c(rel.height.legend, 1)
258		}
259
260	283x	layout_matrix = matrix(c(1, 1, 2:(n_plots + 1)),
261	283x	n_plot_rows + 1, 2, byrow = TRUE)
262	283x	layout(layout_matrix, heights = rel.heights)
263
264	283x	par(mar = c(0.1, 2.1, 0.1, 2.1))
265
266		# Empty plot with legend
267	!	if (!is.null(pred_over)) lty_over <- seq(2, length.out = length(pred_over))
268	283x	else lty_over <- NULL
269	283x	if (pop_curves) {
270	283x	if (is.null(covariates)) {
271	283x	lty_pop <- 1
272	283x	names(lty_pop) <- "Population"
273		} else {
274	!	lty_pop <- 1:nrow(covariates)
275	!	names(lty_pop) <- rownames(covariates)
276		}
277		} else {
278	!	lty_pop <- NULL
279		}
280	283x	n_pop_over <- length(lty_pop) + length(lty_over)
281
282	283x	plot(0, type = "n", axes = FALSE, ann = FALSE)
283	283x	legend("center", bty = "n", ncol = ncol.legend,
284	283x	legend = c(names(lty_pop), names(pred_over), ds_names[i]),
285	283x	lty = c(lty_pop, lty_over, lty_ds),
286	283x	lwd = c(rep(2, n_pop_over), rep(1, length(i))),
287	283x	col = c(rep(1, n_pop_over), col_ds),
288	283x	pch = c(rep(NA, n_pop_over), pch_ds))
289
290	283x	resplot <- match.arg(resplot)
291
292		# Loop plot rows
293	283x	for (plot_row in 1:n_plot_rows) {
294
295	501x	obs_var <- obs_vars[plot_row]
296	501x	observed_row <- subset(observed, name == obs_var)
297
298		# Set ylim to sensible default, or use ymax
299	501x	if (identical(ymax, "auto")) {
300	501x	ylim_row = c(0,
301	501x	max(c(observed_row$value, pred_ds[[obs_var]]), na.rm = TRUE))
302		} else {
303	!	ylim_row = c(0, ymax[plot_row])
304		}
305
306		# Margins for bottom row of plots when we have more than one row
307		# This is the only row that needs to show the x axis legend
308	501x	if (plot_row == n_plot_rows) {
309	283x	par(mar = c(5.1, 4.1, 1.1, 2.1))
310		} else {
311	218x	par(mar = c(3.0, 4.1, 1.1, 2.1))
312		}
313
314	501x	plot(0, type = "n",
315	501x	xlim = xlim, ylim = ylim_row,
316	501x	xlab = xlab, ylab = paste("Residues", obs_var), frame = frame)
317
318	501x	if (!is.null(pred_over)) {
319	!	for (i_over in seq_along(pred_over)) {
320	!	pred_frame <- as.data.frame(pred_over[[i_over]])
321	!	lines(pred_frame$time, pred_frame[[obs_var]],
322	!	lwd = 2, lty = lty_over[i_over])
323		}
324		}
325
326	501x	for (ds_i in seq_along(i)) {
327	4915x	points(subset(observed_row, ds == ds_names[ds_i], c("time", "value")),
328	4915x	col = col_ds[ds_i], pch = pch_ds[ds_i])
329	4915x	lines(subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)),
330	4915x	col = col_ds[ds_i], lty = lty_ds[ds_i])
331		}
332
333	501x	if (pop_curves) {
334	501x	for (cov_i in seq_along(pred_list_pop)) {
335	501x	cov_name <- names(pred_list_pop)[cov_i]
336	501x	lines(
337	501x	pred_list_pop[[cov_i]][, "time"],
338	501x	pred_list_pop[[cov_i]][, obs_var],
339	501x	type = "l", lwd = 2, lty = lty_pop[cov_i])
340		}
341		}
342
343	501x	if (identical(maxabs, "auto")) {
344	283x	maxabs = max(abs(observed_row$residual), na.rm = TRUE)
345		}
346
347	501x	if (identical(resplot, "time")) {
348	!	plot(0, type = "n", xlim = xlim, xlab = "Time",
349	!	ylim = c(-1.2 * maxabs, 1.2 * maxabs),
350	!	ylab = if (standardized) "Standardized residual" else "Residual",
351	!	frame = frame)
352
353	!	abline(h = 0, lty = 2)
354
355	!	for (ds_i in seq_along(i)) {
356	!	points(subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
357	!	col = col_ds[ds_i], pch = pch_ds[ds_i])
358		}
359		}
360
361	501x	if (identical(resplot, "predicted")) {
362	501x	plot(0, type = "n",
363	501x	xlim = c(0, max(pred_ds[[obs_var]])),
364	501x	xlab = "Predicted",
365	501x	ylim = c(-1.2 * maxabs, 1.2 * maxabs),
366	501x	ylab = if (standardized) "Standardized residual" else "Residual",
367	501x	frame = frame)
368
369	501x	abline(h = 0, lty = 2)
370
371	501x	for (ds_i in seq_along(i)) {
372	4915x	observed_row_ds <- merge(
373	4915x	subset(observed_row, ds == ds_names[ds_i], c("time", "residual")),
374	4915x	subset(pred_ds, ds == ds_names[ds_i], c("time", obs_var)))
375	4915x	points(observed_row_ds[c(3, 2)],
376	4915x	col = col_ds[ds_i], pch = pch_ds[ds_i])
377		}
378		}
379		}
380		}

1		#' Function to set up a kinetic model with one or more state variables
2		#'
3		#' This function is usually called using a call to [mkinsub()] for each observed
4		#' variable, specifying the corresponding submodel as well as outgoing pathways
5		#' (see examples).
6		#'
7		#' For the definition of model types and their parameters, the equations given
8		#' in the FOCUS and NAFTA guidance documents are used.
9		#'
10		#' For kinetic models with more than one observed variable, a symbolic solution
11		#' of the system of differential equations is included in the resulting
12		#' mkinmod object in some cases, speeding up the solution.
13		#'
14		#' If a C compiler is found by [pkgbuild::has_compiler()] and there
15		#' is more than one observed variable in the specification, C code is generated
16		#' for evaluating the differential equations, compiled using
17		#' [inline::cfunction()] and added to the resulting mkinmod object.
18		#'
19		#' @param ... For each observed variable, a list as obtained by [mkinsub()]
20		#' has to be specified as an argument (see examples). Currently, single
21		#' first order kinetics "SFO", indeterminate order rate equation kinetics
22		#' "IORE", or single first order with reversible binding "SFORB" are
23		#' implemented for all variables, while "FOMC", "DFOP", "HS" and "logistic"
24		#' can additionally be chosen for the first variable which is assumed to be
25		#' the source compartment.
26		#' Additionally, [mkinsub()] has an argument \code{to}, specifying names of
27		#' variables to which a transfer is to be assumed in the model.
28		#' If the argument \code{use_of_ff} is set to "min"
29		#' and the model for the compartment is "SFO" or "SFORB", an
30		#' additional [mkinsub()] argument can be \code{sink = FALSE}, effectively
31		#' fixing the flux to sink to zero.
32		#' In print.mkinmod, this argument is currently not used.
33		#' @param use_of_ff Specification of the use of formation fractions in the
34		#' model equations and, if applicable, the coefficient matrix. If "max",
35		#' formation fractions are always used (default). If "min", a minimum use of
36		#' formation fractions is made, i.e. each first-order pathway to a metabolite
37		#' has its own rate constant.
38		#' @param speclist The specification of the observed variables and their
39		#' submodel types and pathways can be given as a single list using this
40		#' argument. Default is NULL.