# 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) sum(x$data$residual^2)) names(result$S) <- c("SFO", "IORE", "DFOP") # 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, digits = 3, ...) { cat("Sums of squares:\n") print(x$S) cat("\nCritical sum of squares for checking the SFO model:\n") print(x$S_c) cat("\nParameters:\n") print(x$parameters, digits = digits) t_rep <- .evaluate_nafta_results(x$S, x$S_c, x$distimes, quiet = quiet) cat("\nDTx values:\n") print(signif(x$distimes, digits = digits)) cat("\nRepresentative half-life:\n") print(round(t_rep, 2)) } .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") message("to the terminal degradation rate found with the DFOP model.") message("The representative half-life obtained from the DFOP model may be used") } t_rep <- t_DFOP2 } } return(t_rep) }