Improve interface to experimental version of nlme

The experimental nlme version in my drat repository contains the
variance function structure varConstProp which makes it possible to use
the two-component error model in generalized nonlinear models using
nlme::gnls() and in mixed effects models using nlme::nlme().

9 files changed, 158 insertions, 24 deletions

diff --git a/.travis.yml b/.travis.yml
index b6982328..081a0f83 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -9,8 +9,9 @@ addons:
   apt:
     packages:
       - gcc
-r_packages:
-  - nlme
+repos:
+  CRAN: https://cloud.r-project.org
+  jranke: https://jranke.github.io/drat/
 r_github_packages:
   - jranke/saemixextension
 script:
diff --git a/DESCRIPTION b/DESCRIPTION
index 80ba23d4..b56c77c7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -17,7 +17,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006,
   note that no warranty is implied for correctness of results or fitness for a
   particular purpose.
 Imports: stats, graphics, methods, deSolve, R6, inline, parallel, numDeriv,
-  lmtest, pkgbuild, nlme (>= 3.1-149), purrr, saemix (>= 3.1.9000)
+  lmtest, pkgbuild, nlme (>= 3.1-150.1), purrr, saemix (>= 3.1.9000)
 Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr,
   benchmarkme, tibble, stats4
 License: GPL
diff --git a/R/add_err.R b/R/add_err.R
index d2092a84..8931a281 100644
--- a/R/add_err.R
+++ b/R/add_err.R
@@ -72,7 +72,7 @@
 #' 
 #' @export
 add_err <- function(prediction, sdfunc, secondary = c("M1", "M2"),
-  n = 1000, LOD = 0.1, reps = 2, digits = 1, seed = NA)
+  n = 10, LOD = 0.1, reps = 2, digits = 1, seed = NA)
 {
   if (!is.na(seed)) set.seed(seed)
 
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index a94a26f7..7f7e34e9 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -47,16 +47,16 @@ get_deg_func <- function() {
 #' @examples
 #' ds <- lapply(experimental_data_for_UBA_2019[6:10],
 #'  function(x) subset(x$data[c("name", "time", "value")], name == "parent"))
-#' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1)
+#' f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1)
 #' library(nlme)
-#' endpoints(f[[1]])
-#' f_nlme <- nlme(f)
-#' print(f_nlme)
-#' endpoints(f_nlme)
+#' f_nlme_sfo <- nlme(f["SFO", ])
+#' f_nlme_dfop <- nlme(f["DFOP", ])
+#' AIC(f_nlme_sfo, f_nlme_dfop)
+#' print(f_nlme_dfop)
+#' endpoints(f_nlme_dfop)
 #' \dontrun{
-#'   f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
+#'   f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1))
 #'   update(f_nlme_2, random = parent_0 ~ 1)
-#'   # Test on some real data
 #'   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"),
@@ -100,6 +100,36 @@ get_deg_func <- function() {
 #'
 #'   endpoints(f_nlme_sfo_sfo)
 #'   endpoints(f_nlme_dfop_sfo)
+#'
+#'   if (findFunction("varConstProp")) { # tc error model for nlme available
+#'     # Attempts to fit metabolite kinetics with the tc error model
+#'     #f_2_tc <- 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, quiet = TRUE,
+#'     #   error_model = "tc")
+#'     #f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ], control = list(maxIter = 100))
+#'     #f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ])
+#'     #f_nlme_dfop_sfo_tc <- update(f_nlme_dfop_sfo, weights = varConstProp(),
+#'     #  control = list(sigma = 1, msMaxIter = 100, pnlsMaxIter = 15))
+#'     # Fitting metabolite kinetics with nlme.mmkin and the two-component
+#'     # error model currently does not work, at least not with these data.
+#'
+#'     f_tc <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, error_model = "tc")
+#'     f_nlme_sfo_tc <- nlme(f_tc["SFO", ])
+#'     f_nlme_dfop_tc <- nlme(f_tc["DFOP", ])
+#'     AIC(f_nlme_sfo, f_nlme_sfo_tc, f_nlme_dfop, f_nlme_dfop_tc)
+#'     print(f_nlme_dfop_tc)
+#'   }
+#'   f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo,
+#'    "DFOP-SFO" = m_dfop_sfo),
+#'     ds_2, quiet = TRUE, error_model = "obs")
+#'   f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ])
+#'   # The same with DFOP-SFO does not converge, apparently the variances of
+#'   # parent and A1 are too similar in this case, so that the model is
+#'   # overparameterised
+#'   #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100))
 #' }
 nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
   fixed, random = fixed,
@@ -145,6 +175,24 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
     thisCall[["random"]] <- pdDiag(thisCall[["fixed"]])
   }
 
+  error_model <- model[[1]]$err_mod
+
+  if (missing(weights)) {
+    thisCall[["weights"]] <- switch(error_model,
+      const = NULL,
+      obs = varIdent(form = ~ 1 | name),
+      tc = varConstProp())
+    sigma <- switch(error_model,
+      tc = 1,
+      NULL)
+  }
+
+  control <- thisCall[["control"]]
+  if (error_model == "tc") {
+    control$sigma = 1
+    thisCall[["control"]] <- control
+  }
+
   val <- do.call("nlme.formula", thisCall)
   val$mmkin_orig <- model
   class(val) <- c("nlme.mmkin", "nlme", "lme")
diff --git a/R/sigma_twocomp.R b/R/sigma_twocomp.R
index 1e012d15..e8a92ced 100644
--- a/R/sigma_twocomp.R
+++ b/R/sigma_twocomp.R
@@ -23,7 +23,34 @@
 #' 
 #'   Rocke, David M. and Lorenzato, Stefan (1995) A two-component model for
 #'   measurement error in analytical chemistry. Technometrics 37(2), 176-184.
+#' @examples
+#' times <- c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+#' d_pred <- data.frame(time = times, parent = 100 * exp(- 0.03 * times))
+#' set.seed(123456)
+#' d_syn <- add_err(d_pred, function(y) sigma_twocomp(y, 1, 0.07),
+#'   reps = 2, n = 1)[[1]]
+#' f_nls <- nls(value ~ SSasymp(time, 0, parent_0, lrc), data = d_syn,
+#'  start = list(parent_0 = 100, lrc = -3))
+#' library(nlme)
+#' f_gnls <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
+#'   data = d_syn, na.action = na.omit,
+#'   start = list(parent_0 = 100, lrc = -3))
+#' if (length(findFunction("varConstProp")) > 0) {
+#'   f_gnls_tc <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
+#'     data = d_syn, na.action = na.omit,
+#'     start = list(parent_0 = 100, lrc = -3),
+#'     weights = varConstProp())
+#'   f_gnls_tc_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
+#'     data = d_syn, na.action = na.omit,
+#'     start = list(parent_0 = 100, lrc = -3),
+#'     control = list(sigma = 1),
+#'     weights = varConstProp())
+#' }
+#' f_mkin <- mkinfit("SFO", d_syn, error_model = "const", quiet = TRUE)
+#' f_mkin_tc <- mkinfit("SFO", d_syn, error_model = "tc", quiet = TRUE)
+#' plot_res(f_mkin_tc, standardized = TRUE)
+#' AIC(f_nls, f_gnls, f_gnls_tc, f_gnls_tc_sf, f_mkin, f_mkin_tc)
 #' @export
-  sigma_twocomp <- function(y, sigma_low, rsd_high) {
-    sqrt(sigma_low^2 + y^2 * rsd_high^2)
-  }
+sigma_twocomp <- function(y, sigma_low, rsd_high) {
+  sqrt(sigma_low^2 + y^2 * rsd_high^2)
+}
diff --git a/check.log b/check.log
index 1d00a516..3dfee39f 100644
--- a/check.log
+++ b/check.log
@@ -1,11 +1,11 @@
 * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’
-* using R version 4.0.2 (2020-06-22)
+* using R version 4.0.3 (2020-10-10)
 * using platform: x86_64-pc-linux-gnu (64-bit)
 * using session charset: UTF-8
 * using options ‘--no-tests --as-cran’
 * checking for file ‘mkin/DESCRIPTION’ ... OK
 * checking extension type ... Package
-* this is package ‘mkin’ version ‘0.9.50.3’
+* this is package ‘mkin’ version ‘0.9.50.4’
 * package encoding: UTF-8
 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
 Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
diff --git a/man/add_err.Rd b/man/add_err.Rd
index 9527d508..fe7002ee 100644
--- a/man/add_err.Rd
+++ b/man/add_err.Rd
@@ -8,7 +8,7 @@ add_err(
   prediction,
   sdfunc,
   secondary = c("M1", "M2"),
-  n = 1000,
+  n = 10,
   LOD = 0.1,
   reps = 2,
   digits = 1,
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
index 10c3ec78..0af670a0 100644
--- a/man/nlme.mmkin.Rd
+++ b/man/nlme.mmkin.Rd
@@ -77,16 +77,16 @@ have been obtained by fitting the same model to a list of datasets.
 \examples{
 ds <- lapply(experimental_data_for_UBA_2019[6:10],
  function(x) subset(x$data[c("name", "time", "value")], name == "parent"))
-f <- mmkin("SFO", ds, quiet = TRUE, cores = 1)
+f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1)
 library(nlme)
-endpoints(f[[1]])
-f_nlme <- nlme(f)
-print(f_nlme)
-endpoints(f_nlme)
+f_nlme_sfo <- nlme(f["SFO", ])
+f_nlme_dfop <- nlme(f["DFOP", ])
+AIC(f_nlme_sfo, f_nlme_dfop)
+print(f_nlme_dfop)
+endpoints(f_nlme_dfop)
 \dontrun{
-  f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
+  f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1))
   update(f_nlme_2, random = parent_0 ~ 1)
-  # Test on some real data
   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"),
@@ -130,6 +130,36 @@ endpoints(f_nlme)
 
   endpoints(f_nlme_sfo_sfo)
   endpoints(f_nlme_dfop_sfo)
+
+  if (findFunction("varConstProp")) { # tc error model for nlme available
+    # Attempts to fit metabolite kinetics with the tc error model
+    #f_2_tc <- 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, quiet = TRUE,
+    #   error_model = "tc")
+    #f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ], control = list(maxIter = 100))
+    #f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ])
+    #f_nlme_dfop_sfo_tc <- update(f_nlme_dfop_sfo, weights = varConstProp(),
+    #  control = list(sigma = 1, msMaxIter = 100, pnlsMaxIter = 15))
+    # Fitting metabolite kinetics with nlme.mmkin and the two-component
+    # error model currently does not work, at least not with these data.
+
+    f_tc <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, error_model = "tc")
+    f_nlme_sfo_tc <- nlme(f_tc["SFO", ])
+    f_nlme_dfop_tc <- nlme(f_tc["DFOP", ])
+    AIC(f_nlme_sfo, f_nlme_sfo_tc, f_nlme_dfop, f_nlme_dfop_tc)
+    print(f_nlme_dfop_tc)
+  }
+  f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo,
+   "DFOP-SFO" = m_dfop_sfo),
+    ds_2, quiet = TRUE, error_model = "obs")
+  f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ])
+  # The same with DFOP-SFO does not converge, apparently the variances of
+  # parent and A1 are too similar in this case, so that the model is
+  # overparameterised
+  #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100))
 }
 }
 \seealso{
diff --git a/man/sigma_twocomp.Rd b/man/sigma_twocomp.Rd
index 4e1f7c38..ed79d493 100644
--- a/man/sigma_twocomp.Rd
+++ b/man/sigma_twocomp.Rd
@@ -31,6 +31,34 @@ proposed by Rocke and Lorenzato (1995) can be written in this form as well,
 but assumes approximate lognormal distribution of errors for high values of
 y.
 }
+\examples{
+times <- c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+d_pred <- data.frame(time = times, parent = 100 * exp(- 0.03 * times))
+set.seed(123456)
+d_syn <- add_err(d_pred, function(y) sigma_twocomp(y, 1, 0.07),
+  reps = 2, n = 1)[[1]]
+f_nls <- nls(value ~ SSasymp(time, 0, parent_0, lrc), data = d_syn,
+ start = list(parent_0 = 100, lrc = -3))
+library(nlme)
+f_gnls <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
+  data = d_syn, na.action = na.omit,
+  start = list(parent_0 = 100, lrc = -3))
+if (length(findFunction("varConstProp")) > 0) {
+  f_gnls_tc <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
+    data = d_syn, na.action = na.omit,
+    start = list(parent_0 = 100, lrc = -3),
+    weights = varConstProp())
+  f_gnls_tc_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc),
+    data = d_syn, na.action = na.omit,
+    start = list(parent_0 = 100, lrc = -3),
+    control = list(sigma = 1),
+    weights = varConstProp())
+}
+f_mkin <- mkinfit("SFO", d_syn, error_model = "const", quiet = TRUE)
+f_mkin_tc <- mkinfit("SFO", d_syn, error_model = "tc", quiet = TRUE)
+plot_res(f_mkin_tc, standardized = TRUE)
+AIC(f_nls, f_gnls, f_gnls_tc, f_gnls_tc_sf, f_mkin, f_mkin_tc)
+}
 \references{
 Werner, Mario, Brooks, Samuel H., and Knott, Lancaster B. (1978)
 Additive, Multiplicative, and Mixed Analytical Errors. Clinical Chemistry