Adapt endpoint() to also work for nlme.mmkin objects

8 files changed, 101 insertions, 22 deletions

diff --git a/NEWS.md b/NEWS.md
index a4976488..37800e6d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,6 @@
 # mkin 0.9.49.10 (unreleased)
 
-- 'nlme.mmkin': An nlme method for mmkin row objects and an associated class with plot method
+- 'nlme.mmkin': An nlme method for mmkin row objects and an associated S3 class with print, plot, anova and endpoint methods
 
 - 'mean_degparms, nlme_data, nlme_function': Three new functions to facilitate building nlme models from mmkin row objects
 
diff --git a/R/endpoints.R b/R/endpoints.R
index f7ee483a..586ef9ff 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -7,9 +7,13 @@
 #' are equivalent to the rate constantes of the DFOP model, but with the
 #' advantage that the SFORB model can also be used for metabolites.
 #'
-#' @param fit An object of class \code{\link{mkinfit}}.
+#' @param fit An object of class \code{\link{mkinfit}} or
+#'   \code{\link{nlme.mmkin}}
 #' @importFrom stats optimize
-#' @return A list with the components mentioned above.
+#' @return A list with a matrix of dissipation times named distimes,
+#'   and, if applicable, a vector of formation fractions named ff
+#'   and, if the SFORB model was in use, a vector of eigenvalues
+#'   of these SFORB models, equivalent to DFOP rate constants
 #' @note The function is used internally by \code{\link{summary.mkinfit}}.
 #' @author Johannes Ranke
 #' @keywords manip
@@ -17,6 +21,10 @@
 #'
 #'   fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
 #'   endpoints(fit)
+#'   \dontrun{
+#'     fit_2 <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE)
+#'     endpoints(fit_2)
+#'   }
 #'
 #' @export
 endpoints <- function(fit) {
@@ -25,8 +33,22 @@ endpoints <- function(fit) {
   # 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
   ep <- list()
-  obs_vars <- fit$obs_vars
-  parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
+  if (inherits(fit, "mkinfit")) {
+    mkinmod <- fit$mkinmod
+    parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
+  } else {
+    if (inherits(fit, "nlme.mmkin")) {
+      mkinmod <- fit$mmkin_orig[[1]]$mkinmod
+      bparms.optim <- backtransform_odeparms(fit$coefficients$fixed,
+        mkinmod,
+        transform_rates = fit$mmkin_orig[[1]]$transform_rates,
+        transform_fractions = fit$mmkin_orig[[1]]$transform_fractions)
+      parms.all <- c(bparms.optim, fit$bparms.fixed)
+    } else {
+      stop("Only implemented for mkinfit and nlme.mmkin objects")
+    }
+  }
+  obs_vars <- names(mkinmod$spec)
   ep$ff <- vector()
   ep$SFORB <- vector()
   ep$distimes <- data.frame(
@@ -34,7 +56,7 @@ endpoints <- function(fit) {
     DT90 = rep(NA, length(obs_vars)),
     row.names = obs_vars)
   for (obs_var in obs_vars) {
-    type = names(fit$mkinmod$map[[obs_var]])[1]
+    type = names(mkinmod$map[[obs_var]])[1]
 
     # Get formation fractions if directly fitted, and calculate remaining fraction to sink
     f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE)
@@ -56,7 +78,7 @@ endpoints <- function(fit) {
       k_tot = sum(parms.all[k_names])
       DT50 = log(2)/k_tot
       DT90 = log(10)/k_tot
-      if (fit$mkinmod$use_of_ff == "min") {
+      if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) {
         for (k_name in k_names)
         {
           ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
@@ -78,7 +100,7 @@ endpoints <- function(fit) {
       n = parms.all[paste("N", obs_var, sep = "_")]
       k = k_tot
       # Use the initial concentration of the parent compound
-      source_name = fit$mkinmod$map[[1]][[1]]
+      source_name = mkinmod$map[[1]][[1]]
       c0 = parms.all[paste(source_name, "0", sep = "_")]
       alpha = 1 / (n - 1)
       beta = (c0^(1 - n))/(k * (n - 1))
@@ -86,7 +108,7 @@ endpoints <- function(fit) {
       DT90 = beta * (10^(1/alpha) - 1)
       DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
       ep$distimes[obs_var, c("DT50back")] = DT50_back
-      if (fit$mkinmod$use_of_ff == "min") {
+      if (mkinmod$use_of_ff == "min") {
         for (k_name in k_names)
         {
           ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
@@ -196,6 +218,7 @@ endpoints <- function(fit) {
     }
     ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)
   }
+  if (length(ep$ff) == 0) ep$ff <- NULL
   if (length(ep$SFORB) == 0) ep$SFORB <- NULL
   return(ep)
 }
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index 6e3467ed..1d6c2e75 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -31,8 +31,10 @@
 #'  function(x) subset(x$data[c("name", "time", "value")], name == "parent"))
 #' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1)
 #' library(nlme)
+#' endpoints(f[[1]])
 #' f_nlme <- nlme(f)
 #' print(f_nlme)
+#' endpoints(f_nlme)
 #' f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
 #' update(f_nlme_2, random = parent_0 ~ 1)
 #' \dontrun{
@@ -77,6 +79,9 @@
 #'
 #'   anova(f_nlme_dfop_sfo, f_nlme_fomc_sfo, f_nlme_sfo_sfo)
 #'   anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo) # if we ignore FOMC
+#'
+#'   endpoints(f_nlme_sfo_sfo)
+#'   endpoints(f_nlme_dfop_sfo)
 #' }
 # Code inspired by nlme.nlsList
 nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
diff --git a/docs/news/index.html b/docs/news/index.html
index 5fd97344..cb18664e 100644
--- a/docs/news/index.html
+++ b/docs/news/index.html
@@ -134,7 +134,7 @@
 <a href="#mkin-0-9-49-10-unreleased" class="anchor"></a>mkin 0.9.49.10 (unreleased)<small> Unreleased </small>
 </h1>
 <ul>
-<li><p>‘nlme.mmkin’: An nlme method for mmkin row objects and an associated class with plot method</p></li>
+<li><p>‘nlme.mmkin’: An nlme method for mmkin row objects and an associated S3 class with print, plot, anova and endpoint methods</p></li>
 <li><p>‘mean_degparms, nlme_data, nlme_function’: Three new functions to facilitate building nlme models from mmkin row objects</p></li>
 <li><p>‘endpoints’: Don’t return the SFORB list component if it’s empty. This reduces distraction and complies with the documentation</p></li>
 <li><p>Article in compiled models: Add some platform specific code and suppress warnings about zero values being removed from the FOCUS D dataset</p></li>
diff --git a/docs/reference/endpoints.html b/docs/reference/endpoints.html
index 68af5fcf..ef38c521 100644
--- a/docs/reference/endpoints.html
+++ b/docs/reference/endpoints.html
@@ -154,13 +154,17 @@ advantage that the SFORB model can also be used for metabolites.</p>
     <colgroup><col class="name" /><col class="desc" /></colgroup>
     <tr>
       <th>fit</th>
-      <td><p>An object of class <code><a href='mkinfit.html'>mkinfit</a></code>.</p></td>
+      <td><p>An object of class <code><a href='mkinfit.html'>mkinfit</a></code> or
+<code><a href='nlme.mmkin.html'>nlme.mmkin</a></code></p></td>
     </tr>
     </table>
 
     <h2 class="hasAnchor" id="value"><a class="anchor" href="#value"></a>Value</h2>
 
-    <p>A list with the components mentioned above.</p>
+    <p>A list with a matrix of dissipation times named distimes,
+  and, if applicable, a vector of formation fractions named ff
+  and, if the SFORB model was in use, a vector of eigenvalues
+  of these SFORB models, equivalent to DFOP rate constants</p>
     <h2 class="hasAnchor" id="note"><a class="anchor" href="#note"></a>Note</h2>
 
     <p>The function is used internally by <code><a href='summary.mkinfit.html'>summary.mkinfit</a></code>.</p>
@@ -168,13 +172,24 @@ advantage that the SFORB model can also be used for metabolites.</p>
     <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
     <pre class="examples"><div class='input'>
   <span class='no'>fit</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='st'>"FOMC"</span>, <span class='no'>FOCUS_2006_C</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)
-  <span class='fu'>endpoints</span>(<span class='no'>fit</span>)</div><div class='output co'>#&gt; $ff
-#&gt; logical(0)
-#&gt; 
-#&gt; $distimes
+  <span class='fu'>endpoints</span>(<span class='no'>fit</span>)</div><div class='output co'>#&gt; $distimes
 #&gt;            DT50    DT90 DT50back
 #&gt; parent 1.785233 15.1479 4.559973
-#&gt; </div><div class='input'>
+#&gt; </div><div class='input'>  <span class='co'># \dontrun{</span>
+    <span class='no'>fit_2</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='st'>"SFORB"</span>, <span class='no'>FOCUS_2006_C</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)
+    <span class='fu'>endpoints</span>(<span class='no'>fit_2</span>)</div><div class='output co'>#&gt; $ff
+#&gt; parent_free_sink 
+#&gt;                1 
+#&gt; 
+#&gt; $SFORB
+#&gt; parent_b1 parent_b2 
+#&gt; 0.4595574 0.0178488 
+#&gt; 
+#&gt; $distimes
+#&gt;            DT50     DT90 DT50_parent_b1 DT50_parent_b2
+#&gt; parent 1.886925 21.25106       1.508293       38.83438
+#&gt; </div><div class='input'>  # }
+
 </div></pre>
   </div>
   <div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
diff --git a/docs/reference/nlme.mmkin.html b/docs/reference/nlme.mmkin.html
index 01287dda..9af9cbf7 100644
--- a/docs/reference/nlme.mmkin.html
+++ b/docs/reference/nlme.mmkin.html
@@ -253,7 +253,10 @@ parameters taken from the mmkin object are used</p></td>
  <span class='kw'>function</span>(<span class='no'>x</span>) <span class='fu'><a href='https://rdrr.io/r/base/subset.html'>subset</a></span>(<span class='no'>x</span>$<span class='no'>data</span>[<span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span>(<span class='st'>"name"</span>, <span class='st'>"time"</span>, <span class='st'>"value"</span>)], <span class='no'>name</span> <span class='kw'>==</span> <span class='st'>"parent"</span>))
 <span class='no'>f</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span>(<span class='st'>"SFO"</span>, <span class='no'>ds</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>cores</span> <span class='kw'>=</span> <span class='fl'>1</span>)
 <span class='fu'><a href='https://rdrr.io/r/base/library.html'>library</a></span>(<span class='no'>nlme</span>)
-<span class='no'>f_nlme</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span>(<span class='no'>f</span>)
+<span class='fu'><a href='endpoints.html'>endpoints</a></span>(<span class='no'>f</span><span class='kw'>[[</span><span class='fl'>1</span>]])</div><div class='output co'>#&gt; $distimes
+#&gt;            DT50     DT90
+#&gt; parent 11.96183 39.73634
+#&gt; </div><div class='input'><span class='no'>f_nlme</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span>(<span class='no'>f</span>)
 <span class='fu'><a href='https://rdrr.io/r/base/print.html'>print</a></span>(<span class='no'>f_nlme</span>)</div><div class='output co'>#&gt; Nonlinear mixed-effects model fit by maximum likelihood
 #&gt;   Model: value ~ deg_func(name, time, parent_0, log_k_parent_sink) 
 #&gt;   Data: "Not shown" 
@@ -270,7 +273,10 @@ parameters taken from the mmkin object are used</p></td>
 #&gt; StdDev:  1.30857          1.288591 6.304906
 #&gt; 
 #&gt; Number of Observations: 90
-#&gt; Number of Groups: 5 </div><div class='input'><span class='no'>f_nlme_2</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span>(<span class='no'>f</span>, <span class='kw'>start</span> <span class='kw'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span>(<span class='kw'>parent_0</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>log_k_parent_sink</span> <span class='kw'>=</span> <span class='fl'>0.1</span>))
+#&gt; Number of Groups: 5 </div><div class='input'><span class='fu'><a href='endpoints.html'>endpoints</a></span>(<span class='no'>f_nlme</span>)</div><div class='output co'>#&gt; $distimes
+#&gt;            DT50     DT90
+#&gt; parent 17.51545 58.18505
+#&gt; </div><div class='input'><span class='no'>f_nlme_2</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='https://rdrr.io/pkg/nlme/man/nlme.html'>nlme</a></span>(<span class='no'>f</span>, <span class='kw'>start</span> <span class='kw'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span>(<span class='kw'>parent_0</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>log_k_parent_sink</span> <span class='kw'>=</span> <span class='fl'>0.1</span>))
 <span class='fu'><a href='https://rdrr.io/r/stats/update.html'>update</a></span>(<span class='no'>f_nlme_2</span>, <span class='kw'>random</span> <span class='kw'>=</span> <span class='no'>parent_0</span> ~ <span class='fl'>1</span>)</div><div class='output co'>#&gt; Nonlinear mixed-effects model fit by maximum likelihood
 #&gt;   Model: value ~ deg_func(name, time, parent_0, log_k_parent_sink) 
 #&gt;   Data: "Not shown" 
@@ -410,7 +416,24 @@ parameters taken from the mmkin object are used</p></td>
 #&gt; f_nlme_fomc_sfo     2 11  818.5149  853.0087 -398.2575 1 vs 2  21.33913  &lt;.0001
 #&gt; f_nlme_sfo_sfo      3  9 1085.1821 1113.4042 -533.5910 2 vs 3 270.66712  &lt;.0001</div><div class='input'>  <span class='fu'><a href='https://rdrr.io/r/stats/anova.html'>anova</a></span>(<span class='no'>f_nlme_dfop_sfo</span>, <span class='no'>f_nlme_sfo_sfo</span>) <span class='co'># if we ignore FOMC</span></div><div class='output co'>#&gt;                 Model df       AIC       BIC   logLik   Test L.Ratio p-value
 #&gt; f_nlme_dfop_sfo     1 13  843.8541  884.6194 -408.927                       
-#&gt; f_nlme_sfo_sfo      2  9 1085.1821 1113.4042 -533.591 1 vs 2 249.328  &lt;.0001</div><div class='input'># }
+#&gt; f_nlme_sfo_sfo      2  9 1085.1821 1113.4042 -533.591 1 vs 2 249.328  &lt;.0001</div><div class='input'>
+  <span class='fu'><a href='endpoints.html'>endpoints</a></span>(<span class='no'>f_nlme_sfo_sfo</span>)</div><div class='output co'>#&gt; $ff
+#&gt; parent_sink   parent_A1     A1_sink 
+#&gt;   0.5912435   0.4087565   1.0000000 
+#&gt; 
+#&gt; $distimes
+#&gt;            DT50      DT90
+#&gt; parent 19.13517  63.56565
+#&gt; A1     66.02149 219.31865
+#&gt; </div><div class='input'>  <span class='fu'><a href='endpoints.html'>endpoints</a></span>(<span class='no'>f_nlme_dfop_sfo</span>)</div><div class='output co'>#&gt; $ff
+#&gt;   parent_A1 parent_sink     A1_sink 
+#&gt;   0.2768571   0.7231429   1.0000000 
+#&gt; 
+#&gt; $distimes
+#&gt;             DT50     DT90  DT50_k1 DT50_k2
+#&gt; parent  11.07092 104.6325 4.462389 46.2085
+#&gt; A1     162.30937 539.1801       NA      NA
+#&gt; </div><div class='input'># }
 </div></pre>
   </div>
   <div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
diff --git a/man/endpoints.Rd b/man/endpoints.Rd
index be180737..26b5ee08 100644
--- a/man/endpoints.Rd
+++ b/man/endpoints.Rd
@@ -8,10 +8,14 @@ with mkinfit}
 endpoints(fit)
 }
 \arguments{
-\item{fit}{An object of class \code{\link{mkinfit}}.}
+\item{fit}{An object of class \code{\link{mkinfit}} or
+\code{\link{nlme.mmkin}}}
 }
 \value{
-A list with the components mentioned above.
+A list with a matrix of dissipation times named distimes,
+  and, if applicable, a vector of formation fractions named ff
+  and, if the SFORB model was in use, a vector of eigenvalues
+  of these SFORB models, equivalent to DFOP rate constants
 }
 \description{
 This function calculates DT50 and DT90 values as well as formation fractions
@@ -27,6 +31,10 @@ The function is used internally by \code{\link{summary.mkinfit}}.
 
   fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
   endpoints(fit)
+  \dontrun{
+    fit_2 <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE)
+    endpoints(fit_2)
+  }
 
 }
 \author{
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
index 26dcce66..97d444e4 100644
--- a/man/nlme.mmkin.Rd
+++ b/man/nlme.mmkin.Rd
@@ -79,8 +79,10 @@ 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)
+endpoints(f[[1]])
 f_nlme <- nlme(f)
 print(f_nlme)
+endpoints(f_nlme)
 f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
 update(f_nlme_2, random = parent_0 ~ 1)
 \dontrun{
@@ -125,6 +127,9 @@ update(f_nlme_2, random = parent_0 ~ 1)
 
   anova(f_nlme_dfop_sfo, f_nlme_fomc_sfo, f_nlme_sfo_sfo)
   anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo) # if we ignore FOMC
+
+  endpoints(f_nlme_sfo_sfo)
+  endpoints(f_nlme_dfop_sfo)
 }
 }
 \seealso{