First nlme fits for models with a metabolite

16 files changed, 318 insertions, 82 deletions

diff --git a/man/CAKE_export.Rd b/man/CAKE_export.Rd
index 142b4a75..4bcd8581 100644
--- a/man/CAKE_export.Rd
+++ b/man/CAKE_export.Rd
@@ -4,12 +4,21 @@
 \alias{CAKE_export}
 \title{Export a list of datasets format to a CAKE study file}
 \usage{
-CAKE_export(ds, map = c(parent = "Parent"), links = NA,
-  filename = "CAKE_export.csf", path = ".", overwrite = FALSE,
-  study = "Codlemone aerobic soil degradation", description = "",
-  time_unit = "days", res_unit = "\% AR",
-  comment = "Created using mkin::CAKE_export", date = Sys.Date(),
-  optimiser = "IRLS")
+CAKE_export(
+  ds,
+  map = c(parent = "Parent"),
+  links = NA,
+  filename = "CAKE_export.csf",
+  path = ".",
+  overwrite = FALSE,
+  study = "Codlemone aerobic soil degradation",
+  description = "",
+  time_unit = "days",
+  res_unit = "\% AR",
+  comment = "Created using mkin::CAKE_export",
+  date = Sys.Date(),
+  optimiser = "IRLS"
+)
 }
 \arguments{
 \item{ds}{A named list of datasets in long format as compatible with
diff --git a/man/add_err.Rd b/man/add_err.Rd
index 36b98be9..3452923e 100644
--- a/man/add_err.Rd
+++ b/man/add_err.Rd
@@ -4,8 +4,16 @@
 \alias{add_err}
 \title{Add normally distributed errors to simulated kinetic degradation data}
 \usage{
-add_err(prediction, sdfunc, secondary = c("M1", "M2"), n = 1000,
-  LOD = 0.1, reps = 2, digits = 1, seed = NA)
+add_err(
+  prediction,
+  sdfunc,
+  secondary = c("M1", "M2"),
+  n = 1000,
+  LOD = 0.1,
+  reps = 2,
+  digits = 1,
+  seed = NA
+)
 }
 \arguments{
 \item{prediction}{A prediction from a kinetic model as produced by
diff --git a/man/confint.mkinfit.Rd b/man/confint.mkinfit.Rd
index f2521ccd..81d7a6e8 100644
--- a/man/confint.mkinfit.Rd
+++ b/man/confint.mkinfit.Rd
@@ -4,10 +4,19 @@
 \alias{confint.mkinfit}
 \title{Confidence intervals for parameters of mkinfit objects}
 \usage{
-\method{confint}{mkinfit}(object, parm, level = 0.95, alpha = 1 -
-  level, cutoff, method = c("quadratic", "profile"),
-  transformed = TRUE, backtransform = TRUE,
-  cores = round(detectCores()/2), quiet = FALSE, ...)
+\method{confint}{mkinfit}(
+  object,
+  parm,
+  level = 0.95,
+  alpha = 1 - level,
+  cutoff,
+  method = c("quadratic", "profile"),
+  transformed = TRUE,
+  backtransform = TRUE,
+  cores = round(detectCores()/2),
+  quiet = FALSE,
+  ...
+)
 }
 \arguments{
 \item{object}{An \code{\link{mkinfit}} object}
diff --git a/man/lrtest.mkinfit.Rd b/man/lrtest.mkinfit.Rd
index 84d7bc99..8025b883 100644
--- a/man/lrtest.mkinfit.Rd
+++ b/man/lrtest.mkinfit.Rd
@@ -47,7 +47,7 @@ lrtest(sfo_fit, dfop_fit)
 #lrtest(dfop_fit, error_model = "tc")
 #lrtest(dfop_fit, fixed_parms = c(k2 = 0))
 
-# However, this equivalent syntax works for static help pages
+# However, this equivalent syntax also works for static help pages
 lrtest(dfop_fit, update(dfop_fit, error_model = "tc"))
 lrtest(dfop_fit, update(dfop_fit, fixed_parms = c(k2 = 0)))
 }
diff --git a/man/memkin.Rd b/man/memkin.Rd
new file mode 100644
index 00000000..cda9b468
--- /dev/null
+++ b/man/memkin.Rd
@@ -0,0 +1,72 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/memkin.R
+\name{memkin}
+\alias{memkin}
+\title{Estimation of parameter distributions from mmkin row objects}
+\usage{
+memkin(object, ...)
+}
+\arguments{
+\item{object}{An mmkin row object containing several fits of the same model to different datasets}
+
+\item{...}{Additional arguments passed to \code{\link{nlme}}}
+}
+\value{
+A fitted object of class 'memkin'
+}
+\description{
+This function sets up and attempts to fit a mixed effects model to
+an mmkin row object which is essentially a list of mkinfit objects
+that have been obtained by fitting the same model to a list of
+datasets.
+}
+\examples{
+sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+m_SFO <- mkinmod(parent = mkinsub("SFO"))
+d_SFO_1 <- mkinpredict(m_SFO,
+  c(k_parent_sink = 0.1),
+  c(parent = 98), sampling_times)
+d_SFO_1_long <- mkin_wide_to_long(d_SFO_1, time = "time")
+d_SFO_2 <- mkinpredict(m_SFO,
+  c(k_parent_sink = 0.05),
+  c(parent = 102), sampling_times)
+d_SFO_2_long <- mkin_wide_to_long(d_SFO_2, time = "time")
+d_SFO_3 <- mkinpredict(m_SFO,
+  c(k_parent_sink = 0.02),
+  c(parent = 103), sampling_times)
+d_SFO_3_long <- mkin_wide_to_long(d_SFO_3, time = "time")
+
+d1 <- add_err(d_SFO_1, function(value) 3, n = 1)
+d2 <- add_err(d_SFO_2, function(value) 2, n = 1)
+d3 <- add_err(d_SFO_3, function(value) 4, n = 1)
+ds <- c(d1 = d1, d2 = d2, d3 = d3)
+
+f <- mmkin("SFO", ds)
+x <- memkin(f)
+summary(x)
+
+ds_2 <- lapply(experimental_data_for_UBA_2019[6:10],
+ function(x) x$data[c("name", "time", "value")])
+m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+  A1 = mkinsub("SFO"), use_of_ff = "min")
+m_sfo_sfo_ff <- mkinmod(parent = mkinsub("SFO", "A1"),
+  A1 = mkinsub("SFO"), use_of_ff = "max")
+m_fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+  A1 = mkinsub("SFO"))
+m_dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+  A1 = mkinsub("SFO"))
+
+f_2 <- mmkin(list("SFO-SFO" = m_sfo_sfo,
+ "SFO-SFO-ff" = m_sfo_sfo_ff,
+ "FOMC-SFO" = m_fomc_sfo,
+ "DFOP-SFO" = m_dfop_sfo),
+  ds_2)
+
+f_nlme_sfo_sfo <- memkin(f_2[1, ])
+\dontrun{f_nlme_sfo_sfo_ff <- memkin(f_2[2, ])} # does not converge with maxIter = 50
+f_nlme_fomc_sfo <- memkin(f_2[3, ])
+\dontrun{f_nlme_dfop_sfo <- memkin(f_2[4, ]) # apparently underdetermined}
+anova(f_nlme_sfo_sfo, f_nlme_fomc_sfo)
+# The FOMC variant has a lower AIC and has significantly higher likelihood
+update(f_nlme_fomc_sfo)
+}
diff --git a/man/mkinds.Rd b/man/mkinds.Rd
index 0ea562ed..79eb0167 100644
--- a/man/mkinds.Rd
+++ b/man/mkinds.Rd
@@ -5,9 +5,6 @@
 \alias{mkinds}
 \title{A dataset class for mkin}
 \format{An \code{\link{R6Class}} generator object.}
-\usage{
-mkinds
-}
 \description{
 A dataset class for mkin
 }
@@ -36,3 +33,35 @@ mds <- mkinds$new("FOCUS A", FOCUS_2006_A)
 
 }
 \keyword{datasets}
+\section{Methods}{
+\subsection{Public methods}{
+\itemize{
+\item \href{#method-new}{\code{mkinds$new()}}
+\item \href{#method-clone}{\code{mkinds$clone()}}
+}
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-new"></a>}}
+\subsection{Method \code{new()}}{
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{mkinds$new(title = "", data, time_unit = NA, unit = NA)}\if{html}{\out{</div>}}
+}
+
+}
+\if{html}{\out{<hr>}}
+\if{html}{\out{<a id="method-clone"></a>}}
+\subsection{Method \code{clone()}}{
+The objects of this class are cloneable with this method.
+\subsection{Usage}{
+\if{html}{\out{<div class="r">}}\preformatted{mkinds$clone(deep = FALSE)}\if{html}{\out{</div>}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{<div class="arguments">}}
+\describe{
+\item{\code{deep}}{Whether to make a deep clone.}
+}
+\if{html}{\out{</div>}}
+}
+}
+}
diff --git a/man/mkinerrplot.Rd b/man/mkinerrplot.Rd
index 3c53e7f8..9564ec19 100644
--- a/man/mkinerrplot.Rd
+++ b/man/mkinerrplot.Rd
@@ -4,11 +4,20 @@
 \alias{mkinerrplot}
 \title{Function to plot squared residuals and the error model for an mkin object}
 \usage{
-mkinerrplot(object, obs_vars = names(object$mkinmod$map), xlim = c(0,
-  1.1 * max(object$data$predicted)), xlab = "Predicted",
-  ylab = "Squared residual", maxy = "auto", legend = TRUE,
-  lpos = "topright", col_obs = "auto", pch_obs = "auto",
-  frame = TRUE, ...)
+mkinerrplot(
+  object,
+  obs_vars = names(object$mkinmod$map),
+  xlim = c(0, 1.1 * max(object$data$predicted)),
+  xlab = "Predicted",
+  ylab = "Squared residual",
+  maxy = "auto",
+  legend = TRUE,
+  lpos = "topright",
+  col_obs = "auto",
+  pch_obs = "auto",
+  frame = TRUE,
+  ...
+)
 }
 \arguments{
 \item{object}{A fit represented in an \code{\link{mkinfit}} object.}
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index e58e61e2..45036361 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -8,18 +8,33 @@ Rocke, David M. und Lorenzato, Stefan (1995) A two-component model
   for measurement error in analytical chemistry. Technometrics 37(2), 176-184.
 }
 \usage{
-mkinfit(mkinmod, observed, parms.ini = "auto", state.ini = "auto",
-  err.ini = "auto", fixed_parms = NULL,
-  fixed_initials = names(mkinmod$diffs)[-1], from_max_mean = FALSE,
+mkinfit(
+  mkinmod,
+  observed,
+  parms.ini = "auto",
+  state.ini = "auto",
+  err.ini = "auto",
+  fixed_parms = NULL,
+  fixed_initials = names(mkinmod$diffs)[-1],
+  from_max_mean = FALSE,
   solution_type = c("auto", "analytical", "eigen", "deSolve"),
-  method.ode = "lsoda", use_compiled = "auto",
+  method.ode = "lsoda",
+  use_compiled = "auto",
   control = list(eval.max = 300, iter.max = 200),
-  transform_rates = TRUE, transform_fractions = TRUE, quiet = FALSE,
-  atol = 1e-08, rtol = 1e-10, n.outtimes = 100,
+  transform_rates = TRUE,
+  transform_fractions = TRUE,
+  quiet = FALSE,
+  atol = 1e-08,
+  rtol = 1e-10,
+  n.outtimes = 100,
   error_model = c("const", "obs", "tc"),
-  error_model_algorithm = c("auto", "d_3", "direct", "twostep",
-  "threestep", "fourstep", "IRLS", "OLS"), reweight.tol = 1e-08,
-  reweight.max.iter = 10, trace_parms = FALSE, ...)
+  error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep",
+    "IRLS", "OLS"),
+  reweight.tol = 1e-08,
+  reweight.max.iter = 10,
+  trace_parms = FALSE,
+  ...
+)
 }
 \arguments{
 \item{mkinmod}{A list of class \code{\link{mkinmod}}, containing the kinetic
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index 91f285e2..d2b851b6 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -4,8 +4,13 @@
 \alias{mkinmod}
 \title{Function to set up a kinetic model with one or more state variables}
 \usage{
-mkinmod(..., use_of_ff = "min", speclist = NULL, quiet = FALSE,
-  verbose = FALSE)
+mkinmod(
+  ...,
+  use_of_ff = "min",
+  speclist = NULL,
+  quiet = FALSE,
+  verbose = FALSE
+)
 }
 \arguments{
 \item{...}{For each observed variable, a list has to be specified as an
diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd
index 53f02dea..17d7ef21 100644
--- a/man/mkinpredict.Rd
+++ b/man/mkinpredict.Rd
@@ -6,22 +6,47 @@
 \alias{mkinpredict.mkinfit}
 \title{Produce predictions from a kinetic model using specific parameters}
 \usage{
-mkinpredict(x, odeparms, odeini, outtimes = seq(0, 120, by = 0.1),
-  solution_type = "deSolve", use_compiled = "auto",
-  method.ode = "lsoda", atol = 1e-08, rtol = 1e-10,
-  map_output = TRUE, ...)
-
-\method{mkinpredict}{mkinmod}(x, odeparms = c(k_parent_sink = 0.1),
-  odeini = c(parent = 100), outtimes = seq(0, 120, by = 0.1),
-  solution_type = "deSolve", use_compiled = "auto",
-  method.ode = "lsoda", atol = 1e-08, rtol = 1e-10,
-  map_output = TRUE, ...)
-
-\method{mkinpredict}{mkinfit}(x, odeparms = x$bparms.ode,
-  odeini = x$bparms.state, outtimes = seq(0, 120, by = 0.1),
-  solution_type = "deSolve", use_compiled = "auto",
-  method.ode = "lsoda", atol = 1e-08, rtol = 1e-10,
-  map_output = TRUE, ...)
+mkinpredict(
+  x,
+  odeparms,
+  odeini,
+  outtimes = seq(0, 120, by = 0.1),
+  solution_type = "deSolve",
+  use_compiled = "auto",
+  method.ode = "lsoda",
+  atol = 1e-08,
+  rtol = 1e-10,
+  map_output = TRUE,
+  ...
+)
+
+\method{mkinpredict}{mkinmod}(
+  x,
+  odeparms = c(k_parent_sink = 0.1),
+  odeini = c(parent = 100),
+  outtimes = seq(0, 120, by = 0.1),
+  solution_type = "deSolve",
+  use_compiled = "auto",
+  method.ode = "lsoda",
+  atol = 1e-08,
+  rtol = 1e-10,
+  map_output = TRUE,
+  ...
+)
+
+\method{mkinpredict}{mkinfit}(
+  x,
+  odeparms = x$bparms.ode,
+  odeini = x$bparms.state,
+  outtimes = seq(0, 120, by = 0.1),
+  solution_type = "deSolve",
+  use_compiled = "auto",
+  method.ode = "lsoda",
+  atol = 1e-08,
+  rtol = 1e-10,
+  map_output = TRUE,
+  ...
+)
 }
 \arguments{
 \item{x}{A kinetic model as produced by \code{\link{mkinmod}}, or a kinetic
diff --git a/man/mkinresplot.Rd b/man/mkinresplot.Rd
index 465b3038..2a8b2d41 100644
--- a/man/mkinresplot.Rd
+++ b/man/mkinresplot.Rd
@@ -4,11 +4,21 @@
 \alias{mkinresplot}
 \title{Function to plot residuals stored in an mkin object}
 \usage{
-mkinresplot(object, obs_vars = names(object$mkinmod$map), xlim = c(0,
-  1.1 * max(object$data$time)), standardized = FALSE, xlab = "Time",
+mkinresplot(
+  object,
+  obs_vars = names(object$mkinmod$map),
+  xlim = c(0, 1.1 * max(object$data$time)),
+  standardized = FALSE,
+  xlab = "Time",
   ylab = ifelse(standardized, "Standardized residual", "Residual"),
-  maxabs = "auto", legend = TRUE, lpos = "topright",
-  col_obs = "auto", pch_obs = "auto", frame = TRUE, ...)
+  maxabs = "auto",
+  legend = TRUE,
+  lpos = "topright",
+  col_obs = "auto",
+  pch_obs = "auto",
+  frame = TRUE,
+  ...
+)
 }
 \arguments{
 \item{object}{A fit represented in an \code{\link{mkinfit}} object.}
diff --git a/man/mmkin.Rd b/man/mmkin.Rd
index a763fcdf..4bf07370 100644
--- a/man/mmkin.Rd
+++ b/man/mmkin.Rd
@@ -5,11 +5,16 @@
 \title{Fit one or more kinetic models with one or more state variables to one or
 more datasets}
 \usage{
-mmkin(models = c("SFO", "FOMC", "DFOP"), datasets,
-  cores = round(detectCores()/2), cluster = NULL, ...)
+mmkin(
+  models = c("SFO", "FOMC", "DFOP"),
+  datasets,
+  cores = round(detectCores()/2),
+  cluster = NULL,
+  ...
+)
 }
 \arguments{
-\item{models}{Either a character vector of shorthand names like 
+\item{models}{Either a character vector of shorthand names like
 \code{c("SFO", "FOMC", "DFOP", "HS", "SFORB")}, or an optionally named
 list of \code{\link{mkinmod}} objects.}
 
diff --git a/man/plot.mkinfit.Rd b/man/plot.mkinfit.Rd
index 47cc08a4..c3f3134a 100644
--- a/man/plot.mkinfit.Rd
+++ b/man/plot.mkinfit.Rd
@@ -7,22 +7,47 @@
 \alias{plot_err}
 \title{Plot the observed data and the fitted model of an mkinfit object}
 \usage{
-\method{plot}{mkinfit}(x, fit = x, obs_vars = names(fit$mkinmod$map),
-  xlab = "Time", ylab = "Observed", xlim = range(fit$data$time),
-  ylim = "default", col_obs = 1:length(obs_vars), pch_obs = col_obs,
-  lty_obs = rep(1, length(obs_vars)), add = FALSE, legend = !add,
-  show_residuals = FALSE, show_errplot = FALSE, maxabs = "auto",
-  sep_obs = FALSE, rel.height.middle = 0.9, row_layout = FALSE,
-  lpos = "topright", inset = c(0.05, 0.05), show_errmin = FALSE,
-  errmin_digits = 3, frame = TRUE, ...)
-
-plot_sep(fit, show_errmin = TRUE,
-  show_residuals = ifelse(identical(fit$err_mod, "const"), TRUE,
-  "standardized"), ...)
-
-plot_res(fit, sep_obs = FALSE, show_errmin = sep_obs,
+\method{plot}{mkinfit}(
+  x,
+  fit = x,
+  obs_vars = names(fit$mkinmod$map),
+  xlab = "Time",
+  ylab = "Observed",
+  xlim = range(fit$data$time),
+  ylim = "default",
+  col_obs = 1:length(obs_vars),
+  pch_obs = col_obs,
+  lty_obs = rep(1, length(obs_vars)),
+  add = FALSE,
+  legend = !add,
+  show_residuals = FALSE,
+  show_errplot = FALSE,
+  maxabs = "auto",
+  sep_obs = FALSE,
+  rel.height.middle = 0.9,
+  row_layout = FALSE,
+  lpos = "topright",
+  inset = c(0.05, 0.05),
+  show_errmin = FALSE,
+  errmin_digits = 3,
+  frame = TRUE,
+  ...
+)
+
+plot_sep(
+  fit,
+  show_errmin = TRUE,
+  show_residuals = ifelse(identical(fit$err_mod, "const"), TRUE, "standardized"),
+  ...
+)
+
+plot_res(
+  fit,
+  sep_obs = FALSE,
+  show_errmin = sep_obs,
   standardized = ifelse(identical(fit$err_mod, "const"), FALSE, TRUE),
-  ...)
+  ...
+)
 
 plot_err(fit, sep_obs = FALSE, show_errmin = sep_obs, ...)
 }
diff --git a/man/plot.mmkin.Rd b/man/plot.mmkin.Rd
index 605e458e..f14e0362 100644
--- a/man/plot.mmkin.Rd
+++ b/man/plot.mmkin.Rd
@@ -5,10 +5,19 @@
 \title{Plot model fits (observed and fitted) and the residuals for a row or column
 of an mmkin object}
 \usage{
-\method{plot}{mmkin}(x, main = "auto", legends = 1,
-  resplot = c("time", "errmod"), show_errmin = TRUE,
-  errmin_var = "All data", errmin_digits = 3, cex = 0.7,
-  rel.height.middle = 0.9, ymax = "auto", ...)
+\method{plot}{mmkin}(
+  x,
+  main = "auto",
+  legends = 1,
+  resplot = c("time", "errmod"),
+  show_errmin = TRUE,
+  errmin_var = "All data",
+  errmin_digits = 3,
+  cex = 0.7,
+  rel.height.middle = 0.9,
+  ymax = "auto",
+  ...
+)
 }
 \arguments{
 \item{x}{An object of class \code{\link{mmkin}}, with either one row or one
diff --git a/man/summary.mkinfit.Rd b/man/summary.mkinfit.Rd
index fcc1295f..fabe31d0 100644
--- a/man/summary.mkinfit.Rd
+++ b/man/summary.mkinfit.Rd
@@ -5,11 +5,9 @@
 \alias{print.summary.mkinfit}
 \title{Summary method for class "mkinfit"}
 \usage{
-\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE,
-  alpha = 0.05, ...)
+\method{summary}{mkinfit}(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...)
 
-\method{print}{summary.mkinfit}(x, digits = max(3, getOption("digits") -
-  3), ...)
+\method{print}{summary.mkinfit}(x, digits = max(3, getOption("digits") - 3), ...)
 }
 \arguments{
 \item{object}{an object of class \code{\link{mkinfit}}.}
diff --git a/man/transform_odeparms.Rd b/man/transform_odeparms.Rd
index 5c8c90ba..5257fe12 100644
--- a/man/transform_odeparms.Rd
+++ b/man/transform_odeparms.Rd
@@ -5,11 +5,19 @@
 \alias{backtransform_odeparms}
 \title{Functions to transform and backtransform kinetic parameters for fitting}
 \usage{
-transform_odeparms(parms, mkinmod, transform_rates = TRUE,
-  transform_fractions = TRUE)
+transform_odeparms(
+  parms,
+  mkinmod,
+  transform_rates = TRUE,
+  transform_fractions = TRUE
+)
 
-backtransform_odeparms(transparms, mkinmod, transform_rates = TRUE,
-  transform_fractions = TRUE)
+backtransform_odeparms(
+  transparms,
+  mkinmod,
+  transform_rates = TRUE,
+  transform_fractions = TRUE
+)
 }
 \arguments{
 \item{parms}{Parameters of kinetic models as used in the differential