aboutsummaryrefslogtreecommitdiff
path: root/R/nafta.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/nafta.R')
-rw-r--r--R/nafta.R105
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)
+}

Contact - Imprint