Rudimentary support for setting up nlmixr models

- All degradation models are specified as ODE models. This appears to be
fast enough
- Error models are being translated to nlmixr as close to the mkin error
model as possible. When using the 'saem' backend, it appears not to be
possible to use the same error model for more than one observed variable
- No support yet for models with parallel formation of metabolites, where
the ilr transformation is used in mkin per default
- There is a bug in nlmixr which appears to be triggered if the data are
not balanced, see nlmixrdevelopment/nlmixr#530
- There is a print and a plot method, the summary method is not finished

7 files changed, 334 insertions, 29 deletions

diff --git a/man/mean_degparms.Rd b/man/mean_degparms.Rd
new file mode 100644
index 00000000..92ed4c9d
--- /dev/null
+++ b/man/mean_degparms.Rd
@@ -0,0 +1,27 @@
+% 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{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/nlme.Rd b/man/nlme.Rd
index c367868b..e87b7a00 100644
--- a/man/nlme.Rd
+++ b/man/nlme.Rd
@@ -2,36 +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, test_log_parms = FALSE, conf.level = 0.6)
-
 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?}
-
-\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{
 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{
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
index 2fb0488a..a2b45efa 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,
diff --git a/man/nlmixr.mmkin.Rd b/man/nlmixr.mmkin.Rd
new file mode 100644
index 00000000..86bbdc9f
--- /dev/null
+++ b/man/nlmixr.mmkin.Rd
@@ -0,0 +1,188 @@
+% 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,
+  ...,
+  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",
+  test_log_parms = FALSE,
+  conf.level = 0.6,
+  error_model = object[[1]]$err_mod
+)
+
+nlmixr_data(object, ...)
+}
+\arguments{
+\item{object}{An \link{mmkin} row object containing several fits of the same
+\link{mkinmod} model to different datasets}
+
+\item{est}{Estimation method passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}}
+
+\item{control}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}.}
+
+\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{\dots}{Passed to \link{nlmixr_model}}
+
+\item{x}{An nlmixr.mmkin object to print}
+
+\item{digits}{Number of digits to use for printing}
+
+\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{solution_type}{Possibility to specify the solution type in case the
+automatic choice is not desired}
+}
+\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).
+}
+\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{
+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)
+
+f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem")
+f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei")
+
+f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem")
+f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei")
+
+f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem")
+f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei")
+
+f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem")
+f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei")
+
+f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem")
+f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei")
+
+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", ]))
+
+# nlme is comparable to nlmixr with focei, saem finds a better
+# solution, the two-component error model does not improve it
+plot(f_nlmixr_fomc_saem)
+
+\dontrun{
+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"))
+
+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")
+
+# 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 bcab3e74..d87ca22c 100644
--- a/man/plot.mixed.mmkin.Rd
+++ b/man/plot.mixed.mmkin.Rd
@@ -99,12 +99,17 @@ 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],
diff --git a/man/summary.nlmixr.mmkin.Rd b/man/summary.nlmixr.mmkin.Rd
new file mode 100644
index 00000000..03f0ffb2
--- /dev/null
+++ b/man/summary.nlmixr.mmkin.Rd
@@ -0,0 +1,100 @@
+% 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}
+\title{Summary method for class "nlmixr.mmkin"}
+\usage{
+\method{summary}{nlmixr.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...)
+}
+\arguments{
+\item{object}{an object of class \link{nlmix.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.nlmix.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_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 nlmixr
+f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
+  quiet = TRUE, error_model = "obs", 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")
+#f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei")
+summary(f_nlmixr_dfop_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
index 67cb3cbb..86938d31 100644
--- a/man/summary.saem.mmkin.Rd
+++ b/man/summary.saem.mmkin.Rd
@@ -1,30 +1,32 @@
 % 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}
+% Please edit documentation in R/summary.nlmixr.mmkin.R, R/summary.saem.mmkin.R
+\name{print.summary.saem.mmkin}
 \alias{print.summary.saem.mmkin}
+\alias{summary.saem.mmkin}
 \title{Summary method for class "saem.mmkin"}
 \usage{
+\method{print}{summary.saem.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...)
+
 \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{x}{an object of class \link{summary.saem.mmkin}}
 
-\item{data}{logical, indicating whether the full data should be included in
-the summary.}
+\item{digits}{Number of digits to use for printing}
 
 \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{object}{an object of class \link{saem.mmkin}}
 
-\item{digits}{Number of digits to use for printing}
+\item{data}{logical, indicating whether the full data should be included in
+the summary.}
+
+\item{distimes}{logical, indicating whether DT50 and DT90 values should be
+included.}
 }
 \value{
 The summary function returns a list based on the \link[saemix:SaemixObject-class]{saemix::SaemixObject}