Add a nlme method for mmkin row objects

6 files changed, 270 insertions, 23 deletions

diff --git a/NAMESPACE b/NAMESPACE
index 1b53f101..5f305a3d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,7 @@ S3method(lrtest,mkinfit)
 S3method(lrtest,mmkin)
 S3method(mkinpredict,mkinfit)
 S3method(mkinpredict,mkinmod)
+S3method(nlme,mmkin)
 S3method(nobs,mkinfit)
 S3method(parms,mkinfit)
 S3method(plot,mkinfit)
@@ -91,6 +92,7 @@ importFrom(stats,dist)
 importFrom(stats,dnorm)
 importFrom(stats,lm)
 importFrom(stats,logLik)
+importFrom(stats,na.fail)
 importFrom(stats,nlminb)
 importFrom(stats,nobs)
 importFrom(stats,optimize)
diff --git a/R/nlme.R b/R/nlme.R
index 12a3104c..fafaa7b6 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -1,4 +1,4 @@
-#' Estimation of parameter distributions from mmkin row objects
+#' Helper functions to create nlme models from mmkin row objects
 #'
 #' These functions facilitate setting up a nonlinear mixed effects model for
 #' an mmkin row object. An mmkin row object is essentially a list of mkinfit
@@ -81,14 +81,19 @@
 #'   nlme_f_sfo_sfo <- nlme_function(f_2["SFO-SFO", ])
 #'   nlme_f_sfo_sfo_ff <- nlme_function(f_2["SFO-SFO-ff", ])
 #'   nlme_f_fomc_sfo <- nlme_function(f_2["FOMC-SFO", ])
+#'   assign("nlme_f_sfo_sfo", nlme_f_sfo_sfo, globalenv())
+#'   assign("nlme_f_sfo_sfo_ff", nlme_f_sfo_sfo_ff, globalenv())
+#'   assign("nlme_f_fomc_sfo", nlme_f_fomc_sfo, globalenv())
 #'
-#'   # Allowing for correlations between random effects leads to non-convergence
+#'   # Allowing for correlations between random effects (not shown)
+#'   # leads to non-convergence
 #'   f_nlme_sfo_sfo <- nlme(value ~ nlme_f_sfo_sfo(name, time,
 #'        parent_0, log_k_parent_sink, log_k_parent_A1, log_k_A1_sink),
 #'      data = grouped_data_2,
 #'      fixed = parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1,
 #'      random = pdDiag(parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1),
 #'      start = mean_dp_sfo_sfo)
+#'    # augPred does not see to work on this object, so no plot is shown
 #'
 #'   # The same model fitted with transformed formation fractions does not converge
 #'   f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
@@ -98,14 +103,6 @@
 #'      random = pdDiag(parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1),
 #'      start = mean_dp_sfo_sfo_ff)
 #'
-#'   # It does converge with this version of reduced random effects
-#'   f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
-#'        parent_0, log_k_parent, log_k_A1, f_parent_ilr_1),
-#'      data = grouped_data_2,
-#'      fixed = parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1,
-#'      random = pdDiag(parent_0 + log_k_parent ~ 1),
-#'      start = mean_dp_sfo_sfo_ff)
-#'
 #'   f_nlme_fomc_sfo <- nlme(value ~ nlme_f_fomc_sfo(name, time,
 #'        parent_0, log_alpha, log_beta, log_k_A1, f_parent_ilr_1),
 #'      data = grouped_data_2,
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
new file mode 100644
index 00000000..784ba609
--- /dev/null
+++ b/R/nlme.mmkin.R
@@ -0,0 +1,85 @@
+#' Create an nlme model for an mmkin row object
+#'
+#' This functions sets up a nonlinear mixed effects model for an mmkin row
+#' object. 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.
+#'
+#' @param model An \code{\link{mmkin}} row object.
+#' @param data Ignored, data are taken from the mmkin model
+#' @param fixed Ignored, all degradation parameters fitted in the
+#'   mmkin model are used as fixed parameters
+#' @param random If not specified, all fixed effects are complemented
+#'   with uncorrelated random effects
+#' @param groups See the documentation of nlme
+#' @param start If not specified, mean values of the fitted degradation
+#'   parameters taken from the mmkin object are used
+#' @param correlation See the documentation of nlme
+#' @param weights passed to nlme
+#' @param subset passed to nlme
+#' @param method passed to nlme
+#' @param na.action passed to nlme
+#' @param naPattern passed to nlme
+#' @param control passed to nlme
+#' @param verbose passed to nlme
+#' @importFrom stats na.fail
+#' @return Upon success, a fitted nlme.mmkin object, which is
+#'   an nlme object with additional elements
+#' @export
+#' @seealso \code{\link{nlme_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)
+#' library(nlme)
+#' f_nlme <- nlme(f)
+#' nlme(f, random = parent_0 ~ 1)
+#' f_nlme <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
+#' update(f_nlme, random = parent_0 ~ 1)
+# Code inspired by nlme.nlsList
+nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
+  fixed, random = fixed,
+  groups, start, correlation = NULL, weights = NULL,
+  subset, method = c("ML", "REML"),
+  na.action = na.fail, naPattern,
+  control = list(), verbose= FALSE)
+{
+  if (nrow(model) > 1) stop("Only row objects allowed")
+
+  thisCall <- as.list(match.call())[-1]
+
+  # warn in case of use of arguments that are overriden
+  if (any(!is.na(match(names(thisCall),
+               c("fixed", "data"))))) {
+    warning("'nlme.mmkin' will redefine 'fixed' and 'data'")
+  }
+
+  deg_func <- nlme_function(model)
+  assign("deg_func", deg_func, parent.frame())
+
+  # specify the model formula
+  this_model_text <- paste0("value ~ deg_func(",
+    paste(names(formals(deg_func)), collapse = ", "), ")")
+  this_model <- eval(parse(text = this_model_text))
+  thisCall[["model"]] <- this_model
+
+  mean_dp <- mean_degparms(model)
+  dp_names <- names(mean_dp)
+
+  thisCall[["data"]] <- nlme_data(model)
+
+  if (missing(start)) {
+    thisCall[["start"]] <- mean_dp
+  }
+
+  thisCall[["fixed"]] <- lapply(as.list(dp_names), function(el)
+                                 eval(parse(text = paste(el, 1, sep = "~"))))
+
+  if (missing(random)) {
+    thisCall[["random"]] <- pdDiag(thisCall[["fixed"]])
+  }
+
+  val <- do.call("nlme.formula", thisCall)
+  return(val)
+  class(val) <- c("nlme.mmkin", "nlme", "lme")
+}
+
diff --git a/check.log b/check.log
index 0a75a713..fc06bc0d 100644
--- a/check.log
+++ b/check.log
@@ -5,7 +5,7 @@
 * using options ‘--no-tests --as-cran’
 * checking for file ‘mkin/DESCRIPTION’ ... OK
 * checking extension type ... Package
-* this is package ‘mkin’ version ‘0.9.49.9’
+* this is package ‘mkin’ version ‘0.9.49.10’
 * package encoding: UTF-8
 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
 Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
@@ -41,7 +41,14 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
 * checking S3 generic/method consistency ... OK
 * checking replacement functions ... OK
 * checking foreign function calls ... OK
-* checking R code for possible problems ... OK
+* checking R code for possible problems ... NOTE
+nlme.mmkin: no visible binding for global variable ‘na.fail’
+nlme.mmkin: no visible binding for global variable ‘f’
+Undefined global functions or variables:
+  f na.fail
+Consider adding
+  importFrom("stats", "na.fail")
+to your NAMESPACE file.
 * checking Rd files ... OK
 * checking Rd metadata ... OK
 * checking Rd line widths ... OK
@@ -56,7 +63,91 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
 * checking data for ASCII and uncompressed saves ... OK
 * checking installed files from ‘inst/doc’ ... OK
 * checking files in ‘vignettes’ ... OK
-* checking examples ... OK
+* checking examples ... ERROR
+Running examples in ‘mkin-Ex.R’ failed
+The error most likely occurred in:
+
+> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
+> ### Name: nlme_function
+> ### Title: Estimation of parameter distributions from mmkin row objects
+> ### Aliases: nlme_function mean_degparms nlme_data
+> 
+> ### ** 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, cores = 1, quiet = TRUE)
+> mean_dp <- mean_degparms(f)
+> 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.
+> assign("nlme_f", nlme_f, globalenv())
+> assign("grouped_data", grouped_data, globalenv())
+> 
+> library(nlme)
+> m_nlme <- nlme(value ~ nlme_f(name, time, parent_0, log_k_parent_sink),
++   data = grouped_data,
++   fixed = parent_0 + log_k_parent_sink ~ 1,
++   random = pdDiag(parent_0 + log_k_parent_sink ~ 1),
++   start = mean_dp)
+> summary(m_nlme)
+Nonlinear mixed-effects model fit by maximum likelihood
+  Model: value ~ nlme_f(name, time, parent_0, log_k_parent_sink) 
+ Data: grouped_data 
+       AIC      BIC    logLik
+  253.6377 262.9937 -121.8188
+
+Random effects:
+ Formula: list(parent_0 ~ 1, log_k_parent_sink ~ 1)
+ Level: ds
+ Structure: Diagonal
+        parent_0 log_k_parent_sink Residual
+StdDev: 2.486822         0.6481355 2.368137
+
+Fixed effects: parent_0 + log_k_parent_sink ~ 1 
+                      Value Std.Error DF  t-value p-value
+parent_0          101.43096 1.5960157 44 63.55261       0
+log_k_parent_sink  -3.07614 0.3827294 44 -8.03737       0
+ Correlation: 
+                  prnt_0
+log_k_parent_sink 0.011 
+
+Standardized Within-Group Residuals:
+       Min         Q1        Med         Q3        Max 
+-1.9643731 -0.5430765  0.1659516  0.6160408  1.8082871 
+
+Number of Observations: 48
+Number of Groups: 3 
+> plot(augPred(m_nlme, level = 0:1), layout = c(3, 1))
+> 
+> m_nlme <- nlme(nlme_formula,
++   data = grouped_data,
++   fixed = parent_0 + log_k_parent_sink ~ 1,
++   random = pdDiag(parent_0 + log_k_parent_sink ~ 1),
++   start = mean_dp)
+Error in nlme(nlme_formula, data = grouped_data, fixed = parent_0 + log_k_parent_sink ~  : 
+  object 'nlme_formula' not found
+Execution halted
 * checking for unstated dependencies in ‘tests’ ... OK
 * checking tests ... SKIPPED
 * checking for unstated dependencies in vignettes ... OK
@@ -66,5 +157,8 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: OK
+Status: 1 ERROR, 1 NOTE
+See
+  ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
+for details.
 
diff --git a/man/nlme.Rd b/man/nlme.Rd
index 8e5c2aa0..971ba3f5 100644
--- a/man/nlme.Rd
+++ b/man/nlme.Rd
@@ -101,14 +101,19 @@ plot(augPred(m_nlme, level = 0:1), layout = c(3, 1))
   nlme_f_sfo_sfo <- nlme_function(f_2["SFO-SFO", ])
   nlme_f_sfo_sfo_ff <- nlme_function(f_2["SFO-SFO-ff", ])
   nlme_f_fomc_sfo <- nlme_function(f_2["FOMC-SFO", ])
+  assign("nlme_f_sfo_sfo", nlme_f_sfo_sfo, globalenv())
+  assign("nlme_f_sfo_sfo_ff", nlme_f_sfo_sfo_ff, globalenv())
+  assign("nlme_f_fomc_sfo", nlme_f_fomc_sfo, globalenv())
 
-  # Allowing for correlations between random effects leads to non-convergence
+  # Allowing for correlations between random effects (not shown)
+  # leads to non-convergence
   f_nlme_sfo_sfo <- nlme(value ~ nlme_f_sfo_sfo(name, time,
        parent_0, log_k_parent_sink, log_k_parent_A1, log_k_A1_sink),
      data = grouped_data_2,
      fixed = parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1,
      random = pdDiag(parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1),
      start = mean_dp_sfo_sfo)
+   # augPred does not see to work on this object, so no plot is shown
 
   # The same model fitted with transformed formation fractions does not converge
   f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
@@ -118,14 +123,6 @@ plot(augPred(m_nlme, level = 0:1), layout = c(3, 1))
      random = pdDiag(parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1),
      start = mean_dp_sfo_sfo_ff)
 
-  # It does converge with this version of reduced random effects
-  f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
-       parent_0, log_k_parent, log_k_A1, f_parent_ilr_1),
-     data = grouped_data_2,
-     fixed = parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1,
-     random = pdDiag(parent_0 + log_k_parent ~ 1),
-     start = mean_dp_sfo_sfo_ff)
-
   f_nlme_fomc_sfo <- nlme(value ~ nlme_f_fomc_sfo(name, time,
        parent_0, log_alpha, log_beta, log_k_A1, f_parent_ilr_1),
      data = grouped_data_2,
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
new file mode 100644
index 00000000..5f937488
--- /dev/null
+++ b/man/nlme.mmkin.Rd
@@ -0,0 +1,72 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nlme.mmkin.R
+\name{nlme.mmkin}
+\alias{nlme.mmkin}
+\title{Create an nlme model for an mmkin row object}
+\usage{
+\method{nlme}{mmkin}(
+  model,
+  data = sys.frame(sys.parent()),
+  fixed,
+  random = fixed,
+  groups,
+  start,
+  correlation = NULL,
+  weights = NULL,
+  subset,
+  method = c("ML", "REML"),
+  na.action = na.fail,
+  naPattern,
+  control = list(),
+  verbose = FALSE
+)
+}
+\arguments{
+\item{model}{An \code{\link{mmkin}} row object.}
+
+\item{data}{Ignored, data are taken from the mmkin model}
+
+\item{fixed}{Ignored, all degradation parameters fitted in the
+mmkin model are used as fixed parameters}
+
+\item{random}{If not specified, all fixed effects are complemented
+with uncorrelated random effects}
+
+\item{groups}{See the documentation of nlme}
+
+\item{start}{If not specified, mean values of the fitted degradation
+parameters taken from the mmkin object are used}
+
+\item{correlation}{See the documentation of nlme}
+
+\item{weights}{passed to nlme}
+
+\item{subset}{passed to nlme}
+
+\item{method}{passed to nlme}
+
+\item{na.action}{passed to nlme}
+
+\item{naPattern}{passed to nlme}
+
+\item{control}{passed to nlme}
+
+\item{verbose}{passed to nlme}
+}
+\value{
+Upon success, a fitted nlme.mmkin object, which is
+  an nlme object with additional elements
+}
+\description{
+Create an nlme model for an mmkin row object
+}
+\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)
+library(nlme)
+f_nlme <- nlme(f)
+nlme(f, random = parent_0 ~ 1)
+f_nlme <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
+update(f_nlme, random = parent_0 ~ 1)
+}