1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
#' Evaluate parent kinetics using the NAFTA guidance
#'
#' The function fits the SFO, IORE and DFOP models using \code{\link{mmkin}}
#' and returns an object of class \code{nafta} that has methods for printing
#' and plotting.
#'
#' @param ds A dataframe that must contain one variable called "time" with the
#' time values specified by the \code{time} argument, one column called
#' "name" with the grouping of the observed values, and finally one column of
#' observed values called "value".
#' @param title Optional title of the dataset
#' @param quiet Should the evaluation text be shown?
#' @param \dots Further arguments passed to \code{\link{mmkin}} (not for the
#' printing method).
#' @importFrom stats qf
#' @return An list of class \code{nafta}. The list element named "mmkin" is the
#' \code{\link{mmkin}} object containing the fits of the three models. The
#' list element named "title" contains the title of the dataset used. The
#' list element "data" contains the dataset used in the fits.
#' @author Johannes Ranke
#' @source NAFTA (2011) Guidance for evaluating and calculating degradation
#' kinetics in environmental media. NAFTA Technical Working Group on
#' Pesticides
#' \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/guidance-evaluating-and-calculating-degradation}
#' accessed 2019-02-22
#'
#' US EPA (2015) Standard Operating Procedure for Using the NAFTA Guidance to
#' Calculate Representative Half-life Values and Characterizing Pesticide
#' Degradation
#' \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/standard-operating-procedure-using-nafta-guidance}
#' @examples
#'
#' nafta_evaluation <- nafta(NAFTA_SOP_Appendix_D, cores = 1)
#' print(nafta_evaluation)
#' plot(nafta_evaluation)
#'
#' @export
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 the results of the three models used in the NAFTA scheme.
#'
#' The plots are ordered with increasing complexity of the model in this
#' function (SFO, then IORE, then DFOP).
#'
#' Calls \code{\link{plot.mmkin}}.
#'
#' @param x An object of class \code{\link{nafta}}.
#' @param legend Should a legend be added?
#' @param main Possibility to override the main title of the plot.
#' @param \dots Further arguments passed to \code{\link{plot.mmkin}}.
#' @return The function is called for its side effect.
#' @author Johannes Ranke
#' @export
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 objects
#'
#' Print nafta objects. The results for the three models are printed in the
#' order of increasing model complexity, i.e. SFO, then IORE, and finally DFOP.
#'
#' @param x An \code{\link{nafta}} object.
#' @param digits Number of digits to be used for printing parameters and
#' dissipation times.
#' @rdname nafta
#' @export
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)
}
|
