16 files changed, 810 insertions, 27 deletions
diff --git a/man/dimethenamid_2018.Rd b/man/dimethenamid_2018.Rd
index 31a05a26..0d1265be 100644
--- a/man/dimethenamid_2018.Rd
+++ b/man/dimethenamid_2018.Rd
@@ -11,7 +11,7 @@ An \link{mkindsg} object grouping seven datasets with some meta information
 Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria (2018)
 Renewal Assessment Report Dimethenamid-P Volume 3 - B.8 Environmental fate and behaviour
 Rev. 2 - November 2017
-\href{https://open.efsa.europa.eu}{https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716}
+\url{https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716}
 }
 \usage{
 dimethenamid_2018
@@ -31,5 +31,53 @@ specific pieces of information in the comments.
 }
 \examples{
 print(dimethenamid_2018)
+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
+\dontrun{
+dfop_sfo3_plus <- mkinmod(
+  DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
+  M23 = mkinsub("SFO"),
+  M27 = mkinsub("SFO"),
+  M31 = mkinsub("SFO", "M27", sink = FALSE),
+  quiet = TRUE
+)
+f_dmta_mkin_tc <- mmkin(
+  list("DFOP-SFO3+" = dfop_sfo3_plus),
+  dmta_ds, quiet = TRUE, error_model = "tc")
+nlmixr_model(f_dmta_mkin_tc)
+# The focei fit takes about four minutes on my system
+system.time(
+  f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
+    control = nlmixr::foceiControl(print = 500))
+)
+summary(f_dmta_nlmixr_focei)
+plot(f_dmta_nlmixr_focei)
+# Using saemix takes about 18 minutes
+system.time(
+  f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE)
+)
+
+# nlmixr with est = "saem" is pretty fast with default iteration numbers, most
+# of the time (about 2.5 minutes) is spent for calculating the log likelihood at the end
+# The likelihood calculated for the nlmixr fit is much lower than that found by saemix
+# Also, the trace plot and the plot of the individual predictions is not
+# convincing for the parent. It seems we are fitting an overparameterised
+# model, so the result we get strongly depends on starting parameters and control settings.
+system.time(
+  f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+    control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9))
+)
+traceplot(f_dmta_nlmixr_saem$nm)
+summary(f_dmta_nlmixr_saem)
+plot(f_dmta_nlmixr_saem)
+}
 }
 \keyword{datasets}
diff --git a/man/endpoints.Rd b/man/endpoints.Rd
index 0b225e62..a37ff98d 100644
--- a/man/endpoints.Rd
+++ b/man/endpoints.Rd
@@ -8,8 +8,8 @@ with mkinfit}
 endpoints(fit)
 }
 \arguments{
-\item{fit}{An object of class \link{mkinfit} or \link{nlme.mmkin}
-or another object that has list components
+\item{fit}{An object of class \link{mkinfit}, \link{nlme.mmkin}, \link{saem.mmkin} or
+\link{nlmixr.mmkin}. Or another object that has list components
 mkinmod containing an \link{mkinmod} degradation model, and two numeric vectors,
 bparms.optim and bparms.fixed, that contain parameter values
 for that model.}
@@ -32,8 +32,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
 }
 \note{
-The function is used internally by \link{summary.mkinfit}
-and \link{summary.nlme.mmkin}
+The function is used internally by \link{summary.mkinfit},
+\link{summary.nlme.mmkin} and \link{summary.saem.mmkin}.
 }
 \examples{
 
diff --git a/man/intervals.nlmixr.mmkin.Rd b/man/intervals.nlmixr.mmkin.Rd
new file mode 100644
index 00000000..b1ab5da4
--- /dev/null
+++ b/man/intervals.nlmixr.mmkin.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/intervals.R
+\name{intervals.nlmixr.mmkin}
+\alias{intervals.nlmixr.mmkin}
+\title{Confidence intervals for parameters in nlmixr.mmkin objects}
+\usage{
+\method{intervals}{nlmixr.mmkin}(object, level = 0.95, backtransform = TRUE, ...)
+}
+\arguments{
+\item{object}{The fitted saem.mmkin object}
+
+\item{level}{The confidence level.}
+
+\item{backtransform}{Should we backtransform the parameters where a one to
+one correlation between transformed and backtransformed parameters exists?}
+}
+\value{
+An object with 'intervals.saem.mmkin' and 'intervals.lme' in the
+class attribute
+}
+\description{
+Confidence intervals for parameters in nlmixr.mmkin objects
+}
diff --git a/man/intervals.saem.mmkin.Rd b/man/intervals.saem.mmkin.Rd
new file mode 100644
index 00000000..4a41dd1b
--- /dev/null
+++ b/man/intervals.saem.mmkin.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/intervals.R
+\name{intervals.saem.mmkin}
+\alias{intervals.saem.mmkin}
+\title{Confidence intervals for parameters in saem.mmkin objects}
+\usage{
+\method{intervals}{saem.mmkin}(object, level = 0.95, backtransform = TRUE, ...)
+}
+\arguments{
+\item{object}{The fitted saem.mmkin object}
+
+\item{level}{The confidence level. Must be the default of 0.95 as this is what
+is available in the saemix object}
+
+\item{backtransform}{In case the model was fitted with mkin transformations,
+should we backtransform the parameters where a one to one correlation
+between transformed and backtransformed parameters exists?}
+}
+\value{
+An object with 'intervals.saem.mmkin' and 'intervals.lme' in the
+class attribute
+}
+\description{
+Confidence intervals for parameters in saem.mmkin objects
+}
diff --git a/man/mean_degparms.Rd b/man/mean_degparms.Rd
new file mode 100644
index 00000000..5e2b4b0f
--- /dev/null
+++ b/man/mean_degparms.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mean_degparms.R
+\name{mean_degparms}
+\alias{mean_degparms}
+\title{Calculate mean degradation parameters for an mmkin row object}
+\usage{
+mean_degparms(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
+}
+\arguments{
+\item{object}{An mmkin row object containing several fits of the same model to different datasets}
+
+\item{random}{Should a list with fixed and random effects be returned?}
+
+\item{test_log_parms}{If TRUE, log parameters are only considered in
+the mean calculations if their untransformed counterparts (most likely
+rate constants) pass the t-test for significant difference from zero.}
+
+\item{conf.level}{Possibility to adjust the required confidence level
+for parameter that are tested if requested by 'test_log_parms'.}
+}
+\value{
+If random is FALSE (default), a named vector containing mean values
+of the fitted degradation model parameters. If random is TRUE, a list with
+fixed and random effects, in the format required by the start argument of
+nlme for the case of a single grouping variable ds.
+}
+\description{
+Calculate mean degradation parameters for an mmkin row object
+}
diff --git a/man/mixed.Rd b/man/mixed.Rd
index 8b00382d..95cae364 100644
--- a/man/mixed.Rd
+++ b/man/mixed.Rd
@@ -23,6 +23,10 @@ mixed(object, ...)
 
 \item{digits}{Number of digits to use for printing.}
 }
+\value{
+An object of class 'mixed.mmkin' which has the observed data in a
+single dataframe which is convenient for plotting
+}
 \description{
 Create a mixed effects model from an mmkin row object
 }
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index bf073634..87ce9016 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -118,9 +118,6 @@ variable, specifying the corresponding submodel as well as outgoing pathways
 
 Print mkinmod objects in a way that the user finds his way to get to its
 components.
-
-This is a convenience function to set up the lists used as arguments for
-\code{\link{mkinmod}}.
 }
 \details{
 For the definition of model types and their parameters, the equations given
diff --git a/man/nlme.Rd b/man/nlme.Rd
index 307cca82..e87b7a00 100644
--- a/man/nlme.Rd
+++ b/man/nlme.Rd
@@ -2,29 +2,19 @@
 % Please edit documentation in R/nlme.R
 \name{nlme_function}
 \alias{nlme_function}
-\alias{mean_degparms}
 \alias{nlme_data}
 \title{Helper functions to create nlme models from mmkin row objects}
 \usage{
 nlme_function(object)
 
-mean_degparms(object, random = FALSE)
-
 nlme_data(object)
 }
 \arguments{
 \item{object}{An mmkin row object containing several fits of the same model to different datasets}
-
-\item{random}{Should a list with fixed and random effects be returned?}
 }
 \value{
 A function that can be used with nlme
 
-If random is FALSE (default), a named vector containing mean values
-of the fitted degradation model parameters. If random is TRUE, a list with
-fixed and random effects, in the format required by the start argument of
-nlme for the case of a single grouping variable ds.
-
 A \code{\link{groupedData}} object
 }
 \description{
@@ -60,7 +50,7 @@ grouped_data <- nlme_data(f)
 nlme_f <- nlme_function(f)
 # These assignments are necessary for these objects to be
 # visible to nlme and augPred when evaluation is done by
-# pkgdown to generated the html docs.
+# pkgdown to generate the html docs.
 assign("nlme_f", nlme_f, globalenv())
 assign("grouped_data", grouped_data, globalenv())
 
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
index 2fb0488a..ed58d603 100644
--- a/man/nlme.mmkin.Rd
+++ b/man/nlme.mmkin.Rd
@@ -13,7 +13,7 @@
     paste(el, 1, sep = "~")))),
   random = pdDiag(fixed),
   groups,
-  start = mean_degparms(model, random = TRUE),
+  start = mean_degparms(model, random = TRUE, test_log_parms = TRUE),
   correlation = NULL,
   weights = NULL,
   subset,
@@ -36,10 +36,9 @@
 \item{fixed}{Ignored, all degradation parameters fitted in the
 mmkin model are used as fixed parameters}
 
-\item{random}{If not specified, correlated random effects are set up
-for all optimised degradation model parameters using the log-Cholesky
-parameterization \link[nlme:pdLogChol]{nlme::pdLogChol} that is also the default of
-the generic \link{nlme} method.}
+\item{random}{If not specified, no correlations between random effects are
+set up for the optimised degradation model parameters. This is
+achieved by using the \link[nlme:pdDiag]{nlme::pdDiag} method.}
 
 \item{groups}{See the documentation of nlme}
 
diff --git a/man/nlmixr.mmkin.Rd b/man/nlmixr.mmkin.Rd
new file mode 100644
index 00000000..f9349727
--- /dev/null
+++ b/man/nlmixr.mmkin.Rd
@@ -0,0 +1,222 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nlmixr.R
+\name{nlmixr.mmkin}
+\alias{nlmixr.mmkin}
+\alias{print.nlmixr.mmkin}
+\alias{nlmixr_model}
+\alias{nlmixr_data}
+\title{Fit nonlinear mixed models using nlmixr}
+\usage{
+\method{nlmixr}{mmkin}(
+  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()
+)
+
+\method{print}{nlmixr.mmkin}(x, digits = max(3, getOption("digits") - 3), ...)
+
+nlmixr_model(
+  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
+)
+
+nlmixr_data(object, ...)
+}
+\arguments{
+\item{object}{An \link{mmkin} row object containing several fits of the same
+\link{mkinmod} model to different datasets}
+
+\item{data}{Not used, as the data are extracted from the mmkin row object}
+
+\item{est}{Estimation method passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
+
+\item{control}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
+
+\item{table}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
+
+\item{error_model}{Optional argument to override the error model which is
+being set based on the error model used in the mmkin row object.}
+
+\item{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 \link{mean_degparms}.}
+
+\item{conf.level}{Possibility to adjust the required confidence level
+for parameter that are tested if requested by 'test_log_parms'.}
+
+\item{degparms_start}{Parameter values given as a named numeric vector will
+be used to override the starting values obtained from the 'mmkin' object.}
+
+\item{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.}
+
+\item{\dots}{Passed to \link{nlmixr_model}}
+
+\item{save}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
+
+\item{envir}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
+
+\item{x}{An nlmixr.mmkin object to print}
+
+\item{digits}{Number of digits to use for printing}
+
+\item{add_attributes}{Should the starting values used for degradation model
+parameters and their distribution and for the error model parameters
+be returned as attributes?}
+}
+\value{
+An S3 object of class 'nlmixr.mmkin', containing the fitted
+\link[nlmixr:nlmixr]{nlmixr::nlmixr} object as a list component named 'nm'. The
+object also inherits from 'mixed.mmkin'.
+
+An function defining a model suitable for fitting with \link[nlmixr:nlmixr]{nlmixr::nlmixr}.
+
+An dataframe suitable for use with \link[nlmixr:nlmixr]{nlmixr::nlmixr}
+}
+\description{
+This function uses \code{\link[nlmixr:nlmixr]{nlmixr::nlmixr()}} as a backend for fitting nonlinear mixed
+effects models created from \link{mmkin} row objects using the Stochastic Approximation
+Expectation Maximisation algorithm (SAEM) or First Order Conditional
+Estimation with Interaction (FOCEI).
+}
+\details{
+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 \link{mkinfit}.
+}
+\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)
+}
+}
+\seealso{
+\link{summary.nlmixr.mmkin} \link{plot.mixed.mmkin}
+}
diff --git a/man/plot.mixed.mmkin.Rd b/man/plot.mixed.mmkin.Rd
index 87a82286..d87ca22c 100644
--- a/man/plot.mixed.mmkin.Rd
+++ b/man/plot.mixed.mmkin.Rd
@@ -13,6 +13,8 @@
   xlim = range(x$data$time),
   resplot = c("predicted", "time"),
   pred_over = NULL,
+  test_log_parms = FALSE,
+  conf.level = 0.6,
   ymax = "auto",
   maxabs = "auto",
   ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
@@ -27,7 +29,7 @@
 )
 }
 \arguments{
-\item{x}{An object of class \link{mixed.mmkin}, \link{nlme.mmkin}}
+\item{x}{An object of class \link{mixed.mmkin}, \link{saem.mmkin} or \link{nlme.mmkin}}
 
 \item{i}{A numeric index to select datasets for which to plot the individual predictions,
 in case plots get too large}
@@ -49,6 +51,12 @@ predicted values?}
 \item{pred_over}{Named list of alternative predictions as obtained
 from \link{mkinpredict} with a compatible \link{mkinmod}.}
 
+\item{test_log_parms}{Passed to \link{mean_degparms} in the case of an
+\link{mixed.mmkin} object}
+
+\item{conf.level}{Passed to \link{mean_degparms} in the case of an
+\link{mixed.mmkin} object}
+
 \item{ymax}{Vector of maximum y axis values}
 
 \item{maxabs}{Maximum absolute value of the residuals. This is used for the
@@ -91,9 +99,23 @@ 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))
 }
 }
 \author{
diff --git a/man/reexports.Rd b/man/reexports.Rd
index ccba7567..dfbb76a7 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -1,10 +1,13 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/lrtest.mkinfit.R, R/nlme.mmkin.R
+% Please edit documentation in R/intervals.R, R/lrtest.mkinfit.R,
+%   R/nlme.mmkin.R, R/nlmixr.R
 \docType{import}
 \name{reexports}
 \alias{reexports}
+\alias{intervals}
 \alias{lrtest}
 \alias{nlme}
+\alias{nlmixr}
 \title{Objects exported from other packages}
 \keyword{internal}
 \description{
@@ -14,6 +17,8 @@ below to see their documentation.
 \describe{
   \item{lmtest}{\code{\link[lmtest]{lrtest}}}
 
-  \item{nlme}{\code{\link[nlme]{nlme}}}
+  \item{nlme}{\code{\link[nlme]{intervals}}, \code{\link[nlme]{nlme}}}
+
+  \item{nlmixr}{\code{\link[nlmixr]{nlmixr}}}
 }}
 
diff --git a/man/saem.Rd b/man/saem.Rd
new file mode 100644
index 00000000..00e9aeda
--- /dev/null
+++ b/man/saem.Rd
@@ -0,0 +1,174 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/saem.R
+\name{saem}
+\alias{saem}
+\alias{saem.mmkin}
+\alias{print.saem.mmkin}
+\alias{saemix_model}
+\alias{saemix_data}
+\title{Fit nonlinear mixed models with SAEM}
+\usage{
+saem(object, ...)
+
+\method{saem}{mmkin}(
+  object,
+  transformations = c("mkin", "saemix"),
+  degparms_start = numeric(),
+  test_log_parms = TRUE,
+  conf.level = 0.6,
+  solution_type = "auto",
+  nbiter.saemix = c(300, 100),
+  control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix,
+    save = FALSE, save.graphs = FALSE),
+  fail_with_errors = TRUE,
+  verbose = FALSE,
+  quiet = FALSE,
+  ...
+)
+
+\method{print}{saem.mmkin}(x, digits = max(3, getOption("digits") - 3), ...)
+
+saemix_model(
+  object,
+  solution_type = "auto",
+  transformations = c("mkin", "saemix"),
+  degparms_start = numeric(),
+  test_log_parms = FALSE,
+  verbose = FALSE,
+  ...
+)
+
+saemix_data(object, verbose = FALSE, ...)
+}
+\arguments{
+\item{object}{An \link{mmkin} row object containing several fits of the same
+\link{mkinmod} model to different datasets}
+
+\item{\dots}{Further parameters passed to \link[saemix:saemixModel]{saemix::saemixModel}.}
+
+\item{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. Currently this is only
+supported in cases where the initial concentration of the parent is not fixed,
+SFO or DFOP is used for the parent and there is either no metabolite or one.}
+
+\item{degparms_start}{Parameter values given as a named numeric vector will
+be used to override the starting values obtained from the 'mmkin' object.}
+
+\item{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 \link{mean_degparms}.}
+
+\item{conf.level}{Possibility to adjust the required confidence level
+for parameter that are tested if requested by 'test_log_parms'.}
+
+\item{solution_type}{Possibility to specify the solution type in case the
+automatic choice is not desired}
+
+\item{nbiter.saemix}{Convenience option to increase the number of
+iterations}
+
+\item{control}{Passed to \link[saemix:saemix]{saemix::saemix}.}
+
+\item{fail_with_errors}{Should a failure to compute standard errors
+from the inverse of the Fisher Information Matrix be a failure?}
+
+\item{verbose}{Should we print information about created objects of
+type \link[saemix:SaemixModel-class]{saemix::SaemixModel} and \link[saemix:SaemixData-class]{saemix::SaemixData}?}
+
+\item{quiet}{Should we suppress the messages saemix prints at the beginning
+and the end of the optimisation process?}
+
+\item{x}{An saem.mmkin object to print}
+
+\item{digits}{Number of digits to use for printing}
+}
+\value{
+An S3 object of class 'saem.mmkin', containing the fitted
+\link[saemix:SaemixObject-class]{saemix::SaemixObject} as a list component named 'so'. The
+object also inherits from 'mixed.mmkin'.
+
+An \link[saemix:SaemixModel-class]{saemix::SaemixModel} object.
+
+An \link[saemix:SaemixData-class]{saemix::SaemixData} object.
+}
+\description{
+This function uses \code{\link[saemix:saemix]{saemix::saemix()}} as a backend for fitting nonlinear mixed
+effects models created from \link{mmkin} row objects using the Stochastic Approximation
+Expectation Maximisation algorithm (SAEM).
+}
+\details{
+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 \link{mkinfit}.
+
+Starting values for the fixed effects (population mean parameters, argument
+psi0 of \code{\link[saemix:saemixModel]{saemix::saemixModel()}} are the mean values of the parameters found
+using \link{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", ])
+
+# 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", ])
+compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so)
+
+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",
+#  control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+
+#saemix::compare.saemix(list(
+#  f_saem_dfop_sfo$so,
+#  f_saem_dfop_sfo_deSolve$so))
+
+# 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))
+}
+}
+\seealso{
+\link{summary.saem.mmkin} \link{plot.mixed.mmkin}
+}
diff --git a/man/summary.nlmixr.mmkin.Rd b/man/summary.nlmixr.mmkin.Rd
new file mode 100644
index 00000000..ab8abd5d
--- /dev/null
+++ b/man/summary.nlmixr.mmkin.Rd
@@ -0,0 +1,103 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/summary.nlmixr.mmkin.R
+\name{summary.nlmixr.mmkin}
+\alias{summary.nlmixr.mmkin}
+\alias{print.summary.nlmixr.mmkin}
+\title{Summary method for class "nlmixr.mmkin"}
+\usage{
+\method{summary}{nlmixr.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...)
+
+\method{print}{summary.nlmixr.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...)
+}
+\arguments{
+\item{object}{an object of class \link{nlmixr.mmkin}}
+
+\item{data}{logical, indicating whether the full data should be included in
+the summary.}
+
+\item{verbose}{Should the summary be verbose?}
+
+\item{distimes}{logical, indicating whether DT50 and DT90 values should be
+included.}
+
+\item{\dots}{optional arguments passed to methods like \code{print}.}
+
+\item{x}{an object of class \link{summary.nlmixr.mmkin}}
+
+\item{digits}{Number of digits to use for printing}
+}
+\value{
+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.
+}
+\description{
+Lists model equations, initial parameter values, optimised parameters
+for fixed effects (population), random effects (deviations from the
+population mean) and residual error model, as well as the resulting
+endpoints such as formation fractions and DT50 values. Optionally
+(default is FALSE), the data are listed in full.
+}
+\examples{
+# Generate five datasets following 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)
+}
+
+}
+\author{
+Johannes Ranke for the mkin specific parts
+nlmixr authors for the parts inherited from nlmixr.
+}
diff --git a/man/summary.saem.mmkin.Rd b/man/summary.saem.mmkin.Rd
new file mode 100644
index 00000000..67cb3cbb
--- /dev/null
+++ b/man/summary.saem.mmkin.Rd
@@ -0,0 +1,100 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/summary.saem.mmkin.R
+\name{summary.saem.mmkin}
+\alias{summary.saem.mmkin}
+\alias{print.summary.saem.mmkin}
+\title{Summary method for class "saem.mmkin"}
+\usage{
+\method{summary}{saem.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...)
+
+\method{print}{summary.saem.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...)
+}
+\arguments{
+\item{object}{an object of class \link{saem.mmkin}}
+
+\item{data}{logical, indicating whether the full data should be included in
+the summary.}
+
+\item{verbose}{Should the summary be verbose?}
+
+\item{distimes}{logical, indicating whether DT50 and DT90 values should be
+included.}
+
+\item{\dots}{optional arguments passed to methods like \code{print}.}
+
+\item{x}{an object of class \link{summary.saem.mmkin}}
+
+\item{digits}{Number of digits to use for printing}
+}
+\value{
+The summary function returns a list based on the \link[saemix:SaemixObject-class]{saemix::SaemixObject}
+obtained in the fit, with at least the following additional components
+\item{saemixversion, mkinversion, Rversion}{The saemix, mkin and R versions used}
+\item{date.fit, date.summary}{The dates where the fit and the summary were
+produced}
+\item{diffs}{The differential equations used in the degradation model}
+\item{use_of_ff}{Was maximum or minimum use made of formation fractions}
+\item{data}{The data}
+\item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals}
+\item{confint_back}{Backtransformed parameters, with confidence intervals if available}
+\item{confint_errmod}{Error model parameters with confidence intervals}
+\item{ff}{The estimated formation fractions derived from the fitted
+model.}
+\item{distimes}{The DT50 and DT90 values for each observed variable.}
+\item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
+The print method is called for its side effect, i.e. printing the summary.
+}
+\description{
+Lists model equations, initial parameter values, optimised parameters
+for fixed effects (population), random effects (deviations from the
+population mean) and residual error model, as well as the resulting
+endpoints such as formation fractions and DT50 values. Optionally
+(default is FALSE), the data are listed in full.
+}
+\examples{
+# Generate five datasets following 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 saem
+f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
+  quiet = TRUE, error_model = "tc", cores = 5)
+f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo)
+summary(f_saem_dfop_sfo, data = TRUE)
+}
+
+}
+\author{
+Johannes Ranke for the mkin specific parts
+saemix authors for the parts inherited from saemix.
+}
diff --git a/man/tffm0.Rd b/man/tffm0.Rd
new file mode 100644
index 00000000..46978d5e
--- /dev/null
+++ b/man/tffm0.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/tffm0.R
+\name{tffm0}
+\alias{tffm0}
+\alias{invtffm0}
+\title{Transform formation fractions as in the first published mkin version}
+\usage{
+tffm0(ff)
+
+invtffm0(ff_trans)
+}
+\arguments{
+\item{ff}{Vector of untransformed formation fractions. The sum
+must be smaller or equal to one}
+
+\item{ff_trans}{Vector of transformed formation fractions that can be
+restricted to the interval from 0 to 1}
+}
+\value{
+A vector of the transformed formation fractions
+
+A vector of backtransformed formation fractions for natural use in degradation models
+}
+\description{
+The transformed fractions can be restricted between 0 and 1 in model
+optimisations. Therefore this transformation was used originally in mkin. It
+was later replaced by the \link{ilr} transformation because the ilr transformed
+fractions can assumed to follow normal distribution. As the ilr
+transformation is not available in \link{RxODE} and can therefore not be used in
+the nlmixr modelling language, this transformation is currently used for
+translating mkin models with formation fractions to more than one target
+compartment for fitting with nlmixr in \link{nlmixr_model}. However,
+this implementation cannot be used there, as it is not accessible
+from RxODE.
+}
+\examples{
+ff_example <- c(
+  0.10983681, 0.09035905, 0.08399383
+)
+ff_example_trans <- tffm0(ff_example)
+invtffm0(ff_example_trans)
+}