Remove nlmixr interface for release of mkin 1.1.0

I am postponing my attempts to get the nlmixr interface to CRAN, given
some problems with nlmixr using R-devel under Windows, see
https://github.com/nlmixrdevelopment/nlmixr/issues/596
and
https://github.com/r-hub/rhub/issues/512,
which is fixed by the removal of nlmixr from the testsuite.

For the tests to be more platform independent, the biphasic mixed
effects models test dataset was defined in a way that fitting
should be more robust (less ill-defined).

7 files changed, 5 insertions, 977 deletions

diff --git a/R/endpoints.R b/R/endpoints.R
index 6bf52f07..e81e7a0a 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -10,8 +10,8 @@
 #' 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], [saem.mmkin] or
-#'  [nlmixr.mmkin]. Or another object that has list components
+#' @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.
diff --git a/R/intervals.R b/R/intervals.R
index 8ab2b7ec..258eb4ad 100644
--- a/R/intervals.R
+++ b/R/intervals.R
@@ -95,87 +95,3 @@ intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...
   attr(res, "level") <- level
   return(res)
 }
-
-#' Confidence intervals for parameters in nlmixr.mmkin objects
-#'
-#' @param object The fitted saem.mmkin object
-#' @param level The confidence level.
-#' @param backtransform Should we backtransform the parameters where a one to
-#'   one correlation between transformed and backtransformed parameters exists?
-#' @param \dots For compatibility with the generic method
-#' @importFrom nlme intervals
-#' @return An object with 'intervals.saem.mmkin' and 'intervals.lme' in the
-#'  class attribute
-#' @export
-intervals.nlmixr.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...)
-{
-
-  # Fixed effects
-  mod_vars <- names(object$mkinmod$diffs)
-
-  conf.int <- confint(object$nm)
-  dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start))
-  ndp <- length(dpnames)
-
-  confint_trans <- as.matrix(conf.int[dpnames, c(3, 1, 4)])
-  colnames(confint_trans) <- c("lower", "est.", "upper")
-
-  if (backtransform) {
-    bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
-      object$transform_rates, object$transform_fractions)
-    bpnames <- names(bp)
-
-    # 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 = "_"), dpnames, 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 dpnames) {
-      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
-      }
-    }
-    confint_ret <- confint_back
-  } else {
-    confint_ret <- confint_trans
-  }
-  attr(confint_ret, "label") <- "Fixed effects:"
-
-  # Random effects
-  ranef_ret <- as.matrix(data.frame(lower = NA, 
-      "est." = sqrt(diag(object$nm$omega)), upper = NA))
-  rownames(ranef_ret) <- paste0(gsub("eta\\.", "sd(", rownames(ranef_ret)), ")")
-  attr(ranef_ret, "label") <- "Random effects:"
-
-  # Error model
-  enames <- names(object$nm$sigma)
-  err_ret <- as.matrix(conf.int[enames, c(3, 1, 4)])
-  colnames(err_ret) <- c("lower", "est.", "upper")
-
-  res <- list(
-    fixed = confint_ret,
-    random = ranef_ret,
-    errmod = err_ret
-  )
-  class(res) <- c("intervals.nlmixr.mmkin", "intervals.lme")
-  attr(res, "level") <- level
-  return(res)
-}
diff --git a/R/nlmixr.R b/R/nlmixr.R
deleted file mode 100644
index 5d05f814..00000000
--- a/R/nlmixr.R
+++ /dev/null
@@ -1,584 +0,0 @@
-utils::globalVariables(c("predicted", "std", "ID", "TIME", "CMT", "DV", "IPRED", "IRES", "IWRES"))
-
-#' @export
-nlmixr::nlmixr
-
-#' Fit nonlinear mixed models using nlmixr
-#'
-#' This function uses [nlmixr::nlmixr()] as a backend for fitting nonlinear mixed
-#' effects models created from [mmkin] row objects using the Stochastic Approximation
-#' Expectation Maximisation algorithm (SAEM) or First Order Conditional
-#' Estimation with Interaction (FOCEI).
-#'
-#' 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].
-#'
-#' @importFrom nlmixr nlmixr tableControl
-#' @importFrom dplyr %>%
-#' @param object An [mmkin] row object containing several fits of the same
-#'   [mkinmod] model to different datasets
-#' @param data Not used, the data are extracted from the mmkin row object
-#' @param est Estimation method passed to [nlmixr::nlmixr]
-#' @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 eta_start Standard deviations on the transformed scale 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 data Not used, as the data are extracted from the mmkin row object
-#' @param table Passed to [nlmixr::nlmixr]
-#' @param error_model Optional argument to override the error model which is
-#'   being set based on the error model used in the mmkin row object.
-#' @param control Passed to [nlmixr::nlmixr]
-#' @param \dots Passed to [nlmixr_model]
-#' @param save Passed to [nlmixr::nlmixr]
-#' @param envir Passed to [nlmixr::nlmixr]
-#' @return An S3 object of class 'nlmixr.mmkin', containing the fitted
-#'   [nlmixr::nlmixr] object as a list component named 'nm'. The
-#'   object also inherits from 'mixed.mmkin'.
-#' @seealso [summary.nlmixr.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 <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1)
-#' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc",
-#'   cores = 1, quiet = TRUE)
-#'
-#' library(nlmixr)
-#' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem",
-#'   control = saemControl(print = 0))
-#' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei",
-#'   control = foceiControl(print = 0))
-#'
-#' f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem",
-#'   control = saemControl(print = 0))
-#' f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei",
-#'   control = foceiControl(print = 0))
-#'
-#' f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem",
-#'   control = saemControl(print = 0))
-#' f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei",
-#'   control = foceiControl(print = 0))
-#'
-#' f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem",
-#'   control = saemControl(print = 0))
-#' f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei",
-#'   control = foceiControl(print = 0))
-#'
-#' f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem",
-#'   control = saemControl(print = 0))
-#' f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei",
-#'   control = foceiControl(print = 0))
-#'
-#' AIC(
-#'   f_nlmixr_sfo_saem$nm, f_nlmixr_sfo_focei$nm,
-#'   f_nlmixr_fomc_saem$nm, f_nlmixr_fomc_focei$nm,
-#'   f_nlmixr_dfop_saem$nm, f_nlmixr_dfop_focei$nm,
-#'   f_nlmixr_hs_saem$nm, f_nlmixr_hs_focei$nm,
-#'   f_nlmixr_fomc_saem_tc$nm, f_nlmixr_fomc_focei_tc$nm)
-#'
-#' AIC(nlme(f_mmkin_parent["FOMC", ]))
-#' AIC(nlme(f_mmkin_parent["HS", ]))
-#'
-#' # The FOCEI fit of FOMC with constant variance or the tc error model is best
-#' plot(f_nlmixr_fomc_saem_tc)
-#'
-#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
-#'   A1 = mkinsub("SFO"), quiet = TRUE)
-#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
-#'   A1 = mkinsub("SFO"), quiet = TRUE)
-#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
-#'   A1 = mkinsub("SFO"), quiet = TRUE)
-#'
-#' f_mmkin_const <- mmkin(list(
-#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
-#'   ds, quiet = TRUE, error_model = "const")
-#' f_mmkin_obs <- mmkin(list(
-#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
-#'   ds, quiet = TRUE, error_model = "obs")
-#' f_mmkin_tc <- mmkin(list(
-#'     "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
-#'   ds, quiet = TRUE, error_model = "tc")
-#'
-#' nlmixr_model(f_mmkin_const["SFO-SFO", ])
-#'
-#' # A single constant variance is currently only possible with est = 'focei' in nlmixr
-#' f_nlmixr_sfo_sfo_focei_const <- nlmixr(f_mmkin_const["SFO-SFO", ], est = "focei")
-#' f_nlmixr_fomc_sfo_focei_const <- nlmixr(f_mmkin_const["FOMC-SFO", ], est = "focei")
-#' f_nlmixr_dfop_sfo_focei_const <- nlmixr(f_mmkin_const["DFOP-SFO", ], est = "focei")
-#'
-#' # Variance by variable is supported by 'saem' and 'focei'
-#' f_nlmixr_fomc_sfo_saem_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "saem")
-#' f_nlmixr_fomc_sfo_focei_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "focei")
-#' f_nlmixr_dfop_sfo_saem_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "saem")
-#' f_nlmixr_dfop_sfo_focei_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "focei")
-#'
-#' # Identical two-component error for all variables is only possible with
-#' # est = 'focei' in nlmixr
-#' f_nlmixr_fomc_sfo_focei_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei")
-#' f_nlmixr_dfop_sfo_focei_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei")
-#'
-#' # Two-component error by variable is possible with both estimation methods
-#' # Variance by variable is supported by 'saem' and 'focei'
-#' f_nlmixr_fomc_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "saem",
-#'   error_model = "obs_tc")
-#' f_nlmixr_fomc_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei",
-#'   error_model = "obs_tc")
-#' f_nlmixr_dfop_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "saem",
-#'   error_model = "obs_tc")
-#' f_nlmixr_dfop_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei",
-#'   error_model = "obs_tc")
-#'
-#' AIC(
-#'   f_nlmixr_sfo_sfo_focei_const$nm,
-#'   f_nlmixr_fomc_sfo_focei_const$nm,
-#'   f_nlmixr_dfop_sfo_focei_const$nm,
-#'   f_nlmixr_fomc_sfo_saem_obs$nm,
-#'   f_nlmixr_fomc_sfo_focei_obs$nm,
-#'   f_nlmixr_dfop_sfo_saem_obs$nm,
-#'   f_nlmixr_dfop_sfo_focei_obs$nm,
-#'   f_nlmixr_fomc_sfo_focei_tc$nm,
-#'   f_nlmixr_dfop_sfo_focei_tc$nm,
-#'   f_nlmixr_fomc_sfo_saem_obs_tc$nm,
-#'   f_nlmixr_fomc_sfo_focei_obs_tc$nm,
-#'   f_nlmixr_dfop_sfo_saem_obs_tc$nm,
-#'   f_nlmixr_dfop_sfo_focei_obs_tc$nm
-#' )
-#' # Currently, FOMC-SFO with two-component error by variable fitted by focei gives the
-#' # lowest AIC
-#' plot(f_nlmixr_fomc_sfo_focei_obs_tc)
-#' summary(f_nlmixr_fomc_sfo_focei_obs_tc)
-#'
-#' # Two parallel metabolites
-#' dmta_ds <- lapply(1:7, function(i) {
-#'   ds_i <- dimethenamid_2018$ds[[i]]$data
-#'   ds_i[ds_i$name == "DMTAP", "name"] <-  "DMTA"
-#'   ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
-#'   ds_i
-#' })
-#' names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
-#' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
-#' dmta_ds[["Elliot 1"]] <- NULL
-#' dmta_ds[["Elliot 2"]] <- NULL
-#' sfo_sfo2 <- mkinmod(
-#'   DMTA = mkinsub("SFO", c("M23", "M27")),
-#'   M23 = mkinsub("SFO"),
-#'   M27 = mkinsub("SFO"),
-#'   quiet = TRUE
-#' )
-#' f_dmta_sfo_sfo2 <- mmkin(
-#'   list("SFO-SFO2" = sfo_sfo2),
-#'   dmta_ds, quiet = TRUE, error_model = "obs")
-#' nlmixr_model(f_dmta_sfo_sfo2)
-#' nlmixr_focei_dmta_sfo_sfo2 <- nlmixr(f_dmta_sfo_sfo2, est = "focei")
-#' }
-#' @export
-nlmixr.mmkin <- function(object, data = NULL,
-  est = NULL, control = list(),
-  table = tableControl(),
-  error_model = object[[1]]$err_mod,
-  test_log_parms = TRUE,
-  conf.level = 0.6,
-  degparms_start = "auto",
-  eta_start = "auto",
-  ...,
-  save = NULL,
-  envir = parent.frame()
-)
-{
-  m_nlmixr <- nlmixr_model(object, est = est,
-    error_model = error_model, add_attributes = TRUE,
-    test_log_parms = test_log_parms, conf.level = conf.level,
-    degparms_start = degparms_start, eta_start = eta_start
-  )
-  d_nlmixr <- nlmixr_data(object)
-
-  mean_dp_start <- attr(m_nlmixr, "mean_dp_start")
-  mean_ep_start <- attr(m_nlmixr, "mean_ep_start")
-
-  attributes(m_nlmixr) <- NULL
-
-  fit_time <- system.time({
-    f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est, control = control)
-  })
-
-  if (is.null(f_nlmixr$CMT)) {
-    nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "IPRED", "IRES", "IWRES")])
-    nlmixr_df$CMT <- as.character(object[[1]]$data$variable[1])
-  } else {
-    nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "CMT", "IPRED", "IRES", "IWRES")])
-  }
-
-  return_data <- nlmixr_df %>%
-    dplyr::transmute(ds = ID, name = CMT, time = TIME, value = DV,
-      predicted = IPRED, residual = IRES,
-      std = IRES/IWRES, standardized = IWRES) %>%
-    dplyr::arrange(ds, name, time)
-
-  bparms_optim <- backtransform_odeparms(f_nlmixr$theta,
-    object[[1]]$mkinmod,
-    object[[1]]$transform_rates,
-    object[[1]]$transform_fractions)
-
-  result <- list(
-    mkinmod = object[[1]]$mkinmod,
-    mmkin = object,
-    transform_rates = object[[1]]$transform_rates,
-    transform_fractions = object[[1]]$transform_fractions,
-    nm = f_nlmixr,
-    est = est,
-    time = fit_time,
-    mean_dp_start = mean_dp_start,
-    mean_ep_start = mean_ep_start,
-    bparms.optim = bparms_optim,
-    bparms.fixed = object[[1]]$bparms.fixed,
-    data = return_data,
-    err_mod = error_model,
-    date.fit = date(),
-    nlmixrversion = as.character(utils::packageVersion("nlmixr")),
-    mkinversion = as.character(utils::packageVersion("mkin")),
-    Rversion = paste(R.version$major, R.version$minor, sep=".")
-  )
-
-  class(result) <- c("nlmixr.mmkin", "mixed.mmkin")
-  return(result)
-}
-
-#' @export
-#' @rdname nlmixr.mmkin
-#' @param x An nlmixr.mmkin object to print
-#' @param digits Number of digits to use for printing
-print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
-  cat("Kinetic nonlinear mixed-effects model fit by", x$est, "using nlmixr")
-  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")
-
-  cat("\nLikelihood:\n")
-  print(data.frame(
-      AIC = AIC(x$nm),
-      BIC = BIC(x$nm),
-      logLik = logLik(x$nm),
-      row.names = " "), digits = digits)
-
-  cat("\nFitted parameters:\n")
-  print(x$nm$parFixed, digits = digits)
-
-  invisible(x)
-}
-
-#' @rdname nlmixr.mmkin
-#' @param add_attributes Should the starting values used for degradation model
-#'   parameters and their distribution and for the error model parameters
-#'   be returned as attributes?
-#' @return An function defining a model suitable for fitting with [nlmixr::nlmixr].
-#' @export
-nlmixr_model <- function(object,
-  est = c("saem", "focei"),
-  degparms_start = "auto",
-  eta_start = "auto",
-  test_log_parms = TRUE, conf.level = 0.6,
-  error_model = object[[1]]$err_mod, add_attributes = FALSE)
-{
-  if (nrow(object) > 1) stop("Only row objects allowed")
-  est = match.arg(est)
-
-  mkin_model <- object[[1]]$mkinmod
-  obs_vars <- names(mkin_model$spec)
-
-  if (error_model == object[[1]]$err_mod) {
-
-    if (length(object[[1]]$mkinmod$spec) > 1 & est == "saem") {
-      if (error_model == "const") {
-        message(
-          "Constant variance for more than one variable is not supported for est = 'saem'\n",
-          "Changing the error model to 'obs' (variance by observed variable)")
-        error_model <- "obs"
-      }
-      if (error_model =="tc") {
-        message(
-          "With est = 'saem', a different error model is required for each observed variable",
-          "Changing the error model to 'obs_tc' (Two-component error for each observed variable)")
-        error_model <- "obs_tc"
-      }
-    }
-  }
-
-  degparms_mmkin <- mean_degparms(object,
-    test_log_parms = test_log_parms,
-    conf.level = conf.level, random = TRUE)
-
-  degparms_optim <- degparms_mmkin$fixed
-
-  degparms_optim_ilr_names <- grep("^f_.*_ilr", names(degparms_optim), value = TRUE)
-  obs_vars_ilr <- unique(gsub("f_(.*)_ilr.*$", "\\1", degparms_optim_ilr_names))
-  degparms_optim_noilr <- degparms_optim[setdiff(names(degparms_optim),
-    degparms_optim_ilr_names)]
-
-  degparms_optim_back <- backtransform_odeparms(degparms_optim,
-      object[[1]]$mkinmod,
-      object[[1]]$transform_rates,
-      object[[1]]$transform_fractions)
-
-  if (degparms_start[1] == "auto") {
-    degparms_start <- degparms_optim_noilr
-    for (obs_var_ilr in obs_vars_ilr) {
-      ff_names <- grep(paste0("^f_", obs_var_ilr, "_"),
-        names(degparms_optim_back), value = TRUE)
-      f_tffm0 <- tffm0(degparms_optim_back[ff_names])
-      f_tffm0_qlogis <- qlogis(f_tffm0)
-      names(f_tffm0_qlogis) <- paste0("f_", obs_var_ilr,
-        "_tffm0_", 1:length(f_tffm0), "_qlogis")
-      degparms_start <- c(degparms_start, f_tffm0_qlogis)
-    }
-  }
-
-  if (eta_start[1] == "auto") {
-    eta_start <- degparms_mmkin$eta[setdiff(names(degparms_optim),
-      degparms_optim_ilr_names)]
-    for (obs_var_ilr in obs_vars_ilr) {
-      ff_n <- length(grep(paste0("^f_", obs_var_ilr, "_"),
-        names(degparms_optim_back), value = TRUE))
-      eta_start_ff <- rep(0.3, ff_n)
-      names(eta_start_ff) <- paste0("f_", obs_var_ilr,
-        "_tffm0_", 1:ff_n, "_qlogis")
-      eta_start <- c(eta_start, eta_start_ff)
-    }
-  }
-
-
-  degparms_fixed <- object[[1]]$bparms.fixed
-
-  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)
-
-  # Definition of the model function
-  f <- function(){}
-
-  ini_block <- "ini({"
-
-  # Initial values for all degradation parameters
-  for (parm_name in names(degparms_start)) {
-    # As initials for state variables are not transformed,
-    # we need to modify the name here as we want to
-    # use the original name in the model block
-    ini_block <- paste0(
-      ini_block,
-      parm_name, " = ",
-      as.character(signif(degparms_start[parm_name], 2)),
-      "\n",
-      "eta.", parm_name, " ~ ",
-      as.character(signif(eta_start[parm_name], 2)),
-      "\n"
-    )
-  }
-
-  # Error model parameters
-  error_model_mkin <- object[[1]]$err_mod
-
-  errparm_names_mkin <- names(object[[1]]$errparms)
-  errparms_mkin <- sapply(errparm_names_mkin, function(parm_name) {
-    mean(sapply(object, function(x) x$errparms[parm_name]))
-  })
-
-  sigma_tc_mkin <- errparms_ini <- errparms_mkin[1] +
-    mean(unlist(sapply(object, function(x) x$data$observed)), na.rm = TRUE) *
-      errparms_mkin[2]
-
-  if (error_model == "const") {
-    if (error_model_mkin == "tc") {
-      errparms_ini <- sigma_tc_mkin
-    } else {
-      errparms_ini <- mean(errparms_mkin)
-    }
-    names(errparms_ini) <- "sigma"
-  }
-
-  if (error_model == "obs") {
-    errparms_ini <- switch(error_model_mkin,
-      const = rep(errparms_mkin["sigma"], length(obs_vars)),
-      obs = errparms_mkin,
-      tc = sigma_tc_mkin)
-    names(errparms_ini) <- paste0("sigma_", obs_vars)
-  }
-
-  if (error_model == "tc") {
-    if (error_model_mkin != "tc") {
-      stop("Not supported")
-    } else {
-      errparms_ini <- errparms_mkin
-    }
-  }
-
-  if (error_model == "obs_tc") {
-    if (error_model_mkin != "tc") {
-      stop("Not supported")
-    } else {
-      errparms_ini <- rep(errparms_mkin, length(obs_vars))
-      names(errparms_ini) <- paste0(
-        rep(names(errparms_mkin), length(obs_vars)),
-        "_",
-        rep(obs_vars, each = 2))
-    }
-  }
-
-  for (parm_name in names(errparms_ini)) {
-    ini_block <- paste0(
-      ini_block,
-      parm_name, " = ",
-      as.character(signif(errparms_ini[parm_name], 2)),
-      "\n"
-    )
-  }
-
-  ini_block <- paste0(ini_block, "})")
-
-  body(f)[2] <- parse(text = ini_block)
-
-  model_block <- "model({"
-
-  # Population initial values for the ODE state variables
-  for (parm_name in odeini_optim_parm_names) {
-    model_block <- paste0(
-      model_block,
-      parm_name, "_model = ",
-      parm_name, " + eta.", parm_name, "\n",
-      gsub("(.*)_0", "\\1(0)", parm_name), " = ", parm_name, "_model\n")
-  }
-
-  # Population initial values for log rate constants
-  for (parm_name in grep("^log_", names(degparms_start), value = TRUE)) {
-    model_block <- paste0(
-      model_block,
-      gsub("^log_", "", parm_name), " = ",
-      "exp(", parm_name, " + eta.", parm_name, ")\n")
-  }
-
-  # Population initial values for logit transformed parameters
-  for (parm_name in grep("_qlogis$", names(degparms_start), value = TRUE)) {
-    if (grepl("_tffm0_", parm_name)) {
-      parm_name_new <- gsub("_qlogis$", "", parm_name)
-    } else {
-      parm_name_new <- names(
-      backtransform_odeparms(degparms_start[parm_name],
-        object[[1]]$mkinmod,
-        object[[1]]$transform_rates,
-        object[[1]]$transform_fractions))
-    }
-    model_block <- paste0(
-      model_block,
-      parm_name_new, " = ",
-      "expit(", parm_name, " + eta.", parm_name, ")\n")
-  }
-
-  # Calculate formation fractions from tffm0 transformed values
-  for (obs_var_ilr in obs_vars_ilr) {
-    ff_names <- grep(paste0("^f_", obs_var_ilr, "_"),
-      names(degparms_optim_back), value = TRUE)
-    pattern <- paste0("^f_", obs_var_ilr, "_to_(.*)$")
-    target_vars <- gsub(pattern, "\\1",
-      grep(paste0("^f_", obs_var_ilr, "_to_"), names(degparms_optim_back), value = TRUE))
-    for (i in 1:length(target_vars)) {
-      ff_name <- ff_names[i]
-      ff_line <- paste0(ff_name, " = f_", obs_var_ilr, "_tffm0_", i)
-      if (i > 1) {
-        for (j in (i - 1):1) {
-          ff_line <- paste0(ff_line, " * (1 - f_", obs_var_ilr, "_tffm0_", j , ")")
-        }
-      }
-      model_block <- paste0(
-        model_block,
-        ff_line,
-        "\n"
-      )
-    }
-  }
-
-  # Differential equations
-  model_block <- paste0(
-    model_block,
-    paste(
-      gsub("d_(.*) =", "d/dt(\\1) =", mkin_model$diffs),
-      collapse = "\n"),
-    "\n"
-  )
-
-  # Error model
-  if (error_model == "const") {
-    model_block <- paste0(model_block,
-    paste(paste0(obs_vars, " ~ add(sigma)"), collapse = "\n"))
-  }
-  if (error_model == "obs") {
-    model_block <- paste0(model_block,
-    paste(paste0(obs_vars, " ~ add(sigma_", obs_vars, ")"), collapse = "\n"),
-    "\n")
-  }
-  if (error_model == "tc") {
-    model_block <- paste0(model_block,
-    paste(paste0(obs_vars, " ~ add(sigma_low) + prop(rsd_high)"), collapse = "\n"),
-      "\n")
-  }
-  if (error_model == "obs_tc") {
-    model_block <- paste0(model_block,
-    paste(
-      paste0(obs_vars, " ~ add(sigma_low_", obs_vars, ") + ",
-        "prop(rsd_high_", obs_vars, ")"), collapse = "\n"),
-      "\n")
-  }
-
-  model_block <- paste0(model_block, "})")
-
-  body(f)[3] <- parse(text = model_block)
-
-  if (add_attributes) {
-    attr(f, "mean_dp_start") <- degparms_optim
-    attr(f, "eta_start") <- degparms_mmkin$eta
-    attr(f, "mean_ep_start") <- errparms_ini
-  }
-
-  return(f)
-}
-
-#' @rdname nlmixr.mmkin
-#' @return An dataframe suitable for use with [nlmixr::nlmixr]
-#' @export
-nlmixr_data <- function(object, ...) {
-  if (nrow(object) > 1) stop("Only row objects allowed")
-  d <- lapply(object, function(x) x$data)
-  compartment_map <- 1:length(object[[1]]$mkinmod$spec)
-  names(compartment_map) <- names(object[[1]]$mkinmod$spec)
-  ds_names <- colnames(object)
-
-  ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
-  names(ds_list) <- ds_names
-  ds_nlmixr <- purrr::map_dfr(ds_list, function(x) x, .id = "ds")
-  ds_nlmixr$variable <- as.character(ds_nlmixr$variable)
-  ds_nlmixr_renamed <- data.frame(
-    ID = ds_nlmixr$ds,
-    TIME = ds_nlmixr$time,
-    AMT = 0, EVID = 0,
-    CMT = ds_nlmixr$variable,
-    DV = ds_nlmixr$observed,
-    stringsAsFactors = FALSE)
-
-  return(ds_nlmixr_renamed)
-}
diff --git a/R/saem.R b/R/saem.R
index 26ea1c8d..d3b23861 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -227,7 +227,8 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
   cat("\nFitted parameters:\n")
   conf.int <- x$so@results@conf.int[c("estimate", "lower", "upper")]
   rownames(conf.int) <- x$so@results@conf.int[["name"]]
-  print(conf.int, digits = digits)
+  conf.int.var <- grepl("^Var\\.", rownames(conf.int))
+  print(conf.int[!conf.int.var, ], digits = digits)
 
   invisible(x)
 }
diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R
deleted file mode 100644
index 94d4ed93..00000000
--- a/R/summary.nlmixr.mmkin.R
+++ /dev/null
@@ -1,250 +0,0 @@
-#' Summary method for class "nlmixr.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.
-#'
-#' @importFrom stats confint sd
-#' @param object an object of class [nlmixr.mmkin]
-#' @param x an object of class [summary.nlmixr.mmkin]
-#' @param data logical, indicating whether the full data should be included in
-#'   the summary.
-#' @param verbose Should the summary be verbose?
-#' @param distimes logical, indicating whether DT50 and DT90 values should be
-#'   included.
-#' @param digits Number of digits to use for printing
-#' @param \dots optional arguments passed to methods like \code{print}.
-#' @return The summary function returns a list obtained in the fit, with at
-#'   least the following additional components
-#'   \item{nlmixrversion, mkinversion, Rversion}{The nlmixr, 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_back}{Backtransformed parameters, with confidence intervals if available}
-#'   \item{ff}{The estimated formation fractions derived from the fitted
-#'      model.}
-#'   \item{distimes}{The DT50 and DT90 values for each observed variable.}
-#'   \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
-#'   The print method is called for its side effect, i.e. printing the summary.
-#' @importFrom stats predict vcov
-#' @author Johannes Ranke for the mkin specific parts
-#'   nlmixr authors for the parts inherited from nlmixr.
-#' @examples
-#' # Generate five datasets following DFOP-SFO kinetics
-#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"),
-#'  m1 = mkinsub("SFO"), quiet = TRUE)
-#' set.seed(1234)
-#' k1_in <- rlnorm(5, log(0.1), 0.3)
-#' k2_in <- rlnorm(5, log(0.02), 0.3)
-#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3))
-#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3))
-#' k_m1_in <- rlnorm(5, log(0.02), 0.3)
-#'
-#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) {
-#'   mkinpredict(dfop_sfo,
-#'     c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1),
-#'     c(parent = 100, m1 = 0),
-#'     sampling_times)
-#' }
-#'
-#' ds_mean_dfop_sfo <- lapply(1:5, function(i) {
-#'   mkinpredict(dfop_sfo,
-#'     c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i],
-#'       f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]),
-#'     c(parent = 100, m1 = 0),
-#'     sampling_times)
-#' })
-#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5)
-#'
-#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) {
-#'   add_err(ds,
-#'     sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
-#'     n = 1)[[1]]
-#' })
-#'
-#' \dontrun{
-#' # Evaluate using mmkin and nlmixr
-#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
-#'   quiet = TRUE, error_model = "tc", cores = 5)
-#' f_saemix_dfop_sfo <- mkin::saem(f_mmkin_dfop_sfo)
-#' f_nlme_dfop_sfo <- mkin::nlme(f_mmkin_dfop_sfo)
-#' f_nlmixr_dfop_sfo_saem <- nlmixr(f_mmkin_dfop_sfo, est = "saem")
-#' # The following takes a very long time but gives
-#' f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei")
-#' AIC(f_nlmixr_dfop_sfo_saem$nm, f_nlmixr_dfop_sfo_focei$nm)
-#' summary(f_nlmixr_dfop_sfo_sfo, data = TRUE)
-#' }
-#'
-#' @export
-summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) {
-
-  mod_vars <- names(object$mkinmod$diffs)
-
-  conf.int <- confint(object$nm)
-  dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start))
-  ndp <- length(dpnames)
-
-  confint_trans <- as.matrix(conf.int[dpnames, c(1, 3, 4)])
-  colnames(confint_trans) <- c("est.", "lower", "upper")
-
-  bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
-    object$transform_rates, object$transform_fractions)
-  bpnames <- names(bp)
-
-  # 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 = "_"), dpnames, 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 dpnames) {
-    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
-    }
-  }
-
-  #  Correlation of fixed effects (inspired by summary.nlme)
-  varFix <- vcov(object$nm)
-  stdFix <- sqrt(diag(varFix))
-  object$corFixed <- array(
-    t(varFix/stdFix)/stdFix,
-    dim(varFix),
-    list(dpnames, dpnames))
-
-  object$confint_trans <- confint_trans
-  object$confint_back <- confint_back
-
-  object$date.summary = date()
-  object$use_of_ff = object$mkinmod$use_of_ff
-
-  object$diffs <- object$mkinmod$diffs
-  object$print_data <- data # boolean: Should we print the data?
-
-  names(object$data)[4] <- "observed" # rename value to observed
-
-  object$verbose <- verbose
-
-  object$fixed <- object$mmkin_orig[[1]]$fixed
-  object$AIC = AIC(object$nm)
-  object$BIC = BIC(object$nm)
-  object$logLik = logLik(object$nm)
-
-  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.nlmixr.mmkin")
-  return(object)
-}
-
-#' @rdname summary.nlmixr.mmkin
-#' @export
-print.summary.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
-  cat("nlmixr version used for fitting:   ", x$nlmixrversion, "\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("\nDegradation model predictions using RxODE\n")
-
-  cat("\nFitted in", x$time[["elapsed"]],  "s\n")
-
-  cat("\nVariance model: ")
-  cat(switch(x$err_mod,
-    const = "Constant variance",
-    obs = "Variance unique to each observed variable",
-    tc = "Two-component variance function",
-    obs_tc = "Two-component variance unique to each observed variable"), "\n")
-
-  cat("\nMean of starting values for individual parameters:\n")
-  print(x$mean_dp_start, digits = digits)
-
-  cat("\nMean of starting values for error model parameters:\n")
-  print(x$mean_ep_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")
-  cat("Likelihood calculated by", nlmixr::getOfvType(x$nm), " \n")
-  print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
-      row.names = " "), digits = digits)
-
-  cat("\nOptimised parameters:\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("\nRandom effects (omega):\n")
-  print(x$nm$omega, digits = digits)
-
-  cat("\nVariance model:\n")
-  print(x$nm$sigma, digits = digits)
-
-  cat("\nBacktransformed parameters:\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)
-}
diff --git a/R/tffm0.R b/R/tffm0.R
deleted file mode 100644
index 56283a5d..00000000
--- a/R/tffm0.R
+++ /dev/null
@@ -1,51 +0,0 @@
-#' Transform formation fractions as in the first published mkin version
-#'
-#' This transformation was used originally in mkin, in order to implement a
-#' constraint for the sum of formation fractions to be smaller than 1.  It was
-#' later replaced by the [ilr] transformation because the ilr transformed
-#' fractions can assumed to follow normal distribution. As the ilr
-#' transformation is not available in [RxODE] and can therefore not be used in
-#' the nlmixr modelling language, the original transformation is currently used
-#' for translating mkin models with formation fractions to more than one target
-#' compartment for fitting with nlmixr in [nlmixr_model]. However, this
-#' implementation cannot be used there, as it is not accessible from RxODE.
-#' 
-#' If the transformed formation fractions are restricted to the interval
-#' between 0 and 1, the sum of backtransformed values is restricted
-#' to this interval.
-#'
-#' @param ff Vector of untransformed formation fractions. The sum
-#'   must be smaller or equal to one
-#' @param ff_trans Vector of transformed formation fractions that can be
-#'   restricted to the interval from 0 to 1
-#' @return A vector of the transformed formation fractions
-#' @export
-#' @examples
-#' ff_example <- c(
-#'   0.10983681, 0.09035905, 0.08399383
-#' )
-#' ff_example_trans <- tffm0(ff_example)
-#' invtffm0(ff_example_trans)
-tffm0 <- function(ff) {
-  n <- length(ff)
-  res <- numeric(n)
-  f_remaining <- 1
-  for (i in 1:n) {
-    res[i] <- ff[i]/f_remaining
-    f_remaining <- f_remaining - ff[i]
-  }
-  return(res)
-}
-#' @rdname tffm0
-#' @export
-#' @return A vector of backtransformed formation fractions for natural use in degradation models
-invtffm0 <- function(ff_trans) {
-  n <- length(ff_trans)
-  res <- numeric(n)
-  f_remaining <- 1
-  for (i in 1:n) {
-    res[i] <- ff_trans[i] * f_remaining
-    f_remaining <- f_remaining - res[i]
-  }
-  return(res)
-}
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 174e7c2d..bf988331 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -230,11 +230,7 @@ backtransform_odeparms <- function(transparms, mkinmod,
       if(transform_fractions) {
         if (any(grepl("qlogis", names(trans_f)))) {
           f_tmp  <- plogis(trans_f)
-          if (any(grepl("_tffm0_.*_qlogis$", names(f_tmp)))) {
-            parms[f_names] <- invtffm0(f_tmp)
-          } else {
-            parms[f_names] <- f_tmp
-          }
+          parms[f_names] <- f_tmp
         } else {
           f_tmp <- invilr(trans_f)
           if (spec[[box]]$sink) {