Implement the NAFTA evaluation scheme

1 files changed, 105 insertions, 0 deletions

diff --git a/R/nafta.R b/R/nafta.R
new file mode 100644
index 00000000..102a9ad6
--- /dev/null
+++ b/R/nafta.R
@@ -0,0 +1,105 @@
+# Copyright (C) 2019 Johannes Ranke
+# Contact: jranke@uni-bremen.de
+
+# This file is part of the R package mkin
+
+# mkin is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+nafta <- function(ds, title = NA, quiet = FALSE) {
+  if (length(levels(ds$name)) > 1) {
+    stop("The NAFTA procedure is only defined for decline data for a single compound")
+  }
+  n <- nrow(subset(ds, !is.na(value)))
+  models <- c("SFO", "IORE", "DFOP")
+
+  result <- list(title = title, data = ds)
+  result$mmkin <- mmkin(models, list(ds), quiet = TRUE)
+
+  distimes <- lapply(result$mmkin, function(x) as.numeric(endpoints(x)$distimes["parent", ]))
+
+  result$distimes <- matrix(NA, nrow = 3, ncol = 3,
+    dimnames = list(models, c("DT50", "DT90", "DT50_rep")))
+  result$distimes["SFO", ] <- distimes[[1]][c(1, 2, 1)]
+  result$distimes["IORE", ] <- distimes[[2]][c(1, 2, 3)]
+  result$distimes["DFOP", ] <- distimes[[3]][c(1, 2, 4)]
+
+  # Get parameters with statistics
+  result$parameters <- lapply(result$mmkin, function(x) {
+    summary(x)$bpar[, c(1, 4:6)]
+  })
+  names(result$parameters) <- models
+
+  # Compare the sum of squared residuals (SSR) to the upper bound of the
+  # confidence region of the SSR for the IORE model
+  result$S <- sapply(result$mmkin, function(x) x$ssr)
+  names(result$S) <- models
+  # Equation (3) on p. 3
+  p <- 3
+  result$S["IORE"]
+  result$S_c <- result$S[["IORE"]] * (1 + p/(n - p) * qf(0.5, p, n - p))
+
+  result$t_rep <- .evaluate_nafta_results(result$S, result$S_c, 
+    result$distimes, quiet = quiet)
+
+  class(result) <- "nafta"
+  return(result)
+}
+
+plot.nafta <- function(x, legend = FALSE, main = "auto", ...) {
+  if (main == "auto") {
+    if (is.na(x$title)) main = ""
+    else main = x$title
+  }
+  plot(x$mmkin, ..., legend = legend, main = main)
+}
+
+print.nafta <- function(x, quiet = TRUE, ...) {
+  cat("Parameters:\n")
+  print(x$parameters)
+  t_rep <- .evaluate_nafta_results(x$S, x$S_c, x$distimes, quiet = quiet)
+  cat("\nDTx values:\n")
+  print(round((x$distimes), digits = 0))
+  cat("\nRepresentative half-life:\n")
+  print(t_rep)
+}
+
+.evaluate_nafta_results <- function(S, S_c, distimes, quiet = FALSE) {
+  t_SFO <- distimes["IORE", "DT50"]
+  t_IORE <- distimes["IORE", "DT50_rep"]
+  t_DFOP2 <- distimes["DFOP", "DT50_rep"]
+
+  if (S["SFO"] < S_c) {
+    if (!quiet) {
+      message("S_SFO is lower than the critical value S_c, use the SFO model")
+    }
+    t_rep <- t_SFO
+  } else {
+    if (!quiet) {
+      message("The SFO model is rejected as S_SFO is equal or higher than the critical value S_c")
+    }
+    if (t_IORE < t_DFOP2) {
+      if (!quiet) {
+        message("The half-life obtained from the IORE model may be used")
+      }
+      t_rep <- t_IORE
+    } else {
+      if (!quiet) {
+        message("The representative half-life of the IORE model is longer than the one corresponding to the terminal degradation rate found with the DFOP model.")
+        message("The reprentative half-life obtained from the DFOP model may be used")
+      }
+      t_rep <- t_DFOP2
+    }
+  }
+  return(t_rep)
+}