diff options
| -rw-r--r-- | NAMESPACE | 2 | ||||
| -rw-r--r-- | R/set_nd_nq.R | 164 | ||||
| -rw-r--r-- | log/test.log | 60 | ||||
| -rw-r--r-- | man/set_nd_nq.Rd | 103 | ||||
| -rw-r--r-- | tests/testthat/print_sfo_saem_1.txt | 6 | ||||
| -rw-r--r-- | tests/testthat/summary_hfit_sfo_tc.txt | 2 | ||||
| -rw-r--r-- | tests/testthat/test_set_nd.R | 86 |
7 files changed, 405 insertions, 18 deletions
@@ -129,6 +129,8 @@ export(plot_sep) export(saem) export(saemix_data) export(saemix_model) +export(set_nd_nq) +export(set_nd_nq_focus) export(sigma_twocomp) export(transform_odeparms) import(deSolve) diff --git a/R/set_nd_nq.R b/R/set_nd_nq.R new file mode 100644 index 00000000..37b9a894 --- /dev/null +++ b/R/set_nd_nq.R @@ -0,0 +1,164 @@ +#' Set non-detects and unquantified values in residue series without replicates +#' +#' This function automates replacing unquantified values in residue time and +#' depth series. For time series, the function performs part of the residue +#' processing proposed in the FOCUS kinetics guidance for parent compounds +#' and metabolites. For two-dimensional residue series over time and depth, +#' it automates the proposal of Boesten et al (2015). +#' +#' @param res_raw Character vector of a residue time series, or matrix of +#' residue values with rows representing depth profiles for a specific sampling +#' time, and columns representing time series of residues at the same depth. +#' Values below the limit of detection (lod) have to be coded as "nd", values +#' between the limit of detection and the limit of quantification, if any, have +#' to be coded as "nq". Samples not analysed have to be coded as "na". All +#' values that are not "na", "nd" or "nq" have to be coercible to numeric +#' @param lod Limit of detection (numeric) +#' @param loq Limit of quantification(numeric). Must be specified if the FOCUS rule to +#' stop after the first non-detection is to be applied +#' @param time_zero_presence Do we assume that residues occur at time zero? +#' This only affects samples from the first sampling time that have been +#' reported as "nd" (not detected). +#' @references Boesten, J. J. T. I., van der Linden, A. M. A., Beltman, W. H. +#' J. and Pol, J. W. (2015). Leaching of plant protection products and their +#' transformation products; Proposals for improving the assessment of leaching +#' to groundwater in the Netherlands — Version 2. Alterra report 2630, Alterra +#' Wageningen UR (University & Research centre) +#' @references FOCUS (2014) Generic Guidance for Estimating Persistence and Degradation +#' Kinetics from Environmental Fate Studies on Pesticides in EU Registration, Version 1.1, +#' 18 December 2014, p. 251 +#' @return A numeric vector, if a vector was supplied, or a numeric matrix otherwise +#' @export +#' @examples +#' # FOCUS (2014) p. 75/76 and 131/132 +#' parent_1 <- c(.12, .09, .05, .03, "nd", "nd", "nd", "nd", "nd", "nd") +#' set_nd_nq(parent_1, 0.02) +#' parent_2 <- c(.12, .09, .05, .03, "nd", "nd", .03, "nd", "nd", "nd") +#' set_nd_nq(parent_2, 0.02) +#' set_nd_nq_focus(parent_2, 0.02, loq = 0.05) +#' parent_3 <- c(.12, .09, .05, .03, "nd", "nd", .06, "nd", "nd", "nd") +#' set_nd_nq(parent_3, 0.02) +#' set_nd_nq_focus(parent_3, 0.02, loq = 0.05) +#' metabolite <- c("nd", "nd", "nd", 0.03, 0.06, 0.10, 0.11, 0.10, 0.09, 0.05, 0.03, "nd", "nd") +#' set_nd_nq(metabolite, 0.02) +#' set_nd_nq_focus(metabolite, 0.02, 0.05) +#' # +#' # Boesten et al. (2015), p. 57/58 +#' table_8 <- matrix( +#' c(10, 10, rep("nd", 4), +#' 10, 10, rep("nq", 2), rep("nd", 2), +#' 10, 10, 10, "nq", "nd", "nd", +#' "nq", 10, "nq", rep("nd", 3), +#' "nd", "nq", "nq", rep("nd", 3), +#' rep("nd", 6), rep("nd", 6)), +#' ncol = 6, byrow = TRUE) +#' set_nd_nq(table_8, 0.5, 1.5, time_zero_presence = TRUE) +#' table_10 <- matrix( +#' c(10, 10, rep("nd", 4), +#' 10, 10, rep("nd", 4), +#' 10, 10, 10, rep("nd", 3), +#' "nd", 10, rep("nd", 4), +#' rep("nd", 18)), +#' ncol = 6, byrow = TRUE) +#' set_nd_nq(table_10, 0.5, time_zero_presence = TRUE) +set_nd_nq <- function(res_raw, lod, loq = NA, time_zero_presence = FALSE) { + if (!is.character(res_raw)) { + stop("Please supply a vector or a matrix of character values") + } + if (is.vector(res_raw)) { + was_vector <- TRUE + res_raw <- as.matrix(res_raw) + } else { + was_vector <- FALSE + if (!is.matrix(res_raw)) { + stop("Please supply a vector or a matrix of character values") + } + } + nq <- 0.5 * (loq + lod) + nda <- 0.5 * lod # not detected but adjacent to detection + res_raw[res_raw == "nq"] <- nq + + if (!time_zero_presence) { + for (j in 1:ncol(res_raw)) { + if (res_raw[1, j] == "nd") res_raw[1, j] <- "na" + } + } + res_raw[res_raw == "na"] <- NA + + not_nd_na <- function(value) !(grepl("nd", value) | is.na(value)) + + for (i in 1:nrow(res_raw)) { + for (j in 1:ncol(res_raw)) { + if (!is.na(res_raw[i, j]) && res_raw[i, j] == "nd") { + if (i > 1) { # check earlier sample in same layer + if (not_nd_na(res_raw[i - 1, j])) res_raw[i, j] <- "nda" + } + if (i < nrow(res_raw)) { # check later sample + if (not_nd_na(res_raw[i + 1, j])) res_raw[i, j] <- "nda" + } + if (j > 1) { # check above sample at the same time + if (not_nd_na(res_raw[i, j - 1])) res_raw[i, j] <- "nda" + } + if (j < ncol(res_raw)) { # check sample below at the same time + if (not_nd_na(res_raw[i, j + 1])) res_raw[i, j] <- "nda" + } + } + } + } + res_raw[res_raw == "nda"] <- nda + res_raw[res_raw == "nd"] <- NA + + result <- as.numeric(res_raw) + dim(result) <- dim(res_raw) + dimnames(result) <- dimnames(res_raw) + if (was_vector) result <- as.vector(result) + return(result) +} + +#' @describeIn set_nd_nq Set non-detects in residue time series according to FOCUS rules +#' @param set_first_sample_nd Should the first sample be set to "first_sample_nd_value" +#' in case it is a non-detection? +#' @param first_sample_nd_value Value to be used for the first sample if it is a non-detection +#' @param ignore_below_loq_after_first_nd Should we ignore values below the LOQ after the first +#' non-detection that occurs after the quantified values? +#' @export +set_nd_nq_focus <- function(res_raw, lod, loq = NA, + set_first_sample_nd = TRUE, first_sample_nd_value = 0, + ignore_below_loq_after_first_nd = TRUE) +{ + + if (!is.vector(res_raw)) stop("FOCUS rules are only specified for one-dimensional time series") + + if (ignore_below_loq_after_first_nd & is.na(loq)) { + stop("You need to specify an LOQ") + } + + n <- length(res_raw) + if (ignore_below_loq_after_first_nd) { + for (i in 3:n) { + if (!res_raw[i - 2] %in% c("na", "nd")) { + if (res_raw[i - 1] == "nd") { + res_remaining <- res_raw[i:n] + res_remaining_unquantified <- ifelse(res_remaining == "na", TRUE, + ifelse(res_remaining == "nd", TRUE, + ifelse(res_remaining == "nq", TRUE, + ifelse(suppressWarnings(as.numeric(res_remaining)) < loq, TRUE, FALSE)))) + res_remaining_numeric <- suppressWarnings(as.numeric(res_remaining)) + res_remaining_below_loq <- ifelse(res_remaining == "nq", TRUE, + ifelse(!is.na(res_remaining_numeric) & res_remaining_numeric < loq, TRUE, FALSE)) + if (all(res_remaining_unquantified)) { + res_raw[i:n] <- ifelse(res_remaining_below_loq, "nd", res_remaining) + } + } + } + } + } + + result <- set_nd_nq(res_raw, lod = lod, loq = loq) + + if (set_first_sample_nd) { + if (res_raw[1] == "nd") result[1] <- first_sample_nd_value + } + + return(result) +} diff --git a/log/test.log b/log/test.log index efa30d17..73fbfcf2 100644 --- a/log/test.log +++ b/log/test.log @@ -1,24 +1,55 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.4s] +✔ | 5 | Analytical solutions for coupled models [3.3s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE ✔ | 12 | Confidence intervals and p-values [1.0s] -✔ | 1 12 | Dimethenamid data from 2018 [31.9s] +✔ | 1 12 | Dimethenamid data from 2018 [32.2s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data Reason: Fitting this ODE model with saemix takes about 15 minutes on my system ──────────────────────────────────────────────────────────────────────────────── -✔ | 14 | Error model fitting [4.9s] +✔ | 14 | Error model fitting [5.1s] ✔ | 5 | Time step normalisation ✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s] ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s] ✔ | 4 | Test fitting the decline of metabolites from their maximum [0.4s] ✔ | 1 | Fitting the logistic model [0.2s] -✔ | 7 | Batch fitting and diagnosing hierarchical kinetic models [14.5s] -✔ | 1 12 | Nonlinear mixed-effects models [0.3s] +✖ | 1 6 | Batch fitting and diagnosing hierarchical kinetic models [14.5s] ──────────────────────────────────────────────────────────────────────────────── +Failure (test_mhmkin.R:41:3): Multiple hierarchical kinetic models can be fitted and diagnosed +Results have changed from known value recorded in 'summary_hfit_sfo_tc.txt'. + + old | new +[45] "" | "" [45] +[46] "Variance model:" | "Variance model:" [46] +[47] " est. lower upper" | " est. lower upper" [47] +[48] "a.1 0.91 0.64 1.17" - "a.1 0.90 0.64 1.17" [48] +[49] "b.1 0.05 0.04 0.06" | "b.1 0.05 0.04 0.06" [49] +[50] "" | "" [50] +[51] "Backtransformed parameters:" | "Backtransformed parameters:" [51] +──────────────────────────────────────────────────────────────────────────────── +✖ | 1 1 11 | Nonlinear mixed-effects models [0.3s] +──────────────────────────────────────────────────────────────────────────────── +Failure (test_mixed.R:8:3): Print methods work +Results have changed from known value recorded in 'print_sfo_saem_1.txt'. + +old[7:19] vs new[7:19] + "" + "Likelihood computed by importance sampling" + " AIC BIC logLik" +- " 1311 1315 -649" ++ " 1312 1316 -650" + "" + "Fitted parameters:" + " estimate lower upper" +- "parent_0 1e+02 99.13 1e+02" ++ "parent_0 1e+02 99.28 1e+02" + "k_parent 4e-02 0.03 4e-02" + "a.1 9e-01 0.75 1e+00" +and 3 more ... + Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve ──────────────────────────────────────────────────────────────────────────────── @@ -26,25 +57,26 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 10 | Special cases of mkinfit calls [0.4s] ✔ | 3 | mkinfit features [0.7s] ✔ | 8 | mkinmod model generation and printing [0.2s] -✔ | 3 | Model predictions with mkinpredict [0.3s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.8s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.3s] -✔ | 16 | Plotting [10.0s] +✔ | 3 | Model predictions with mkinpredict [0.4s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.9s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.6s] +✔ | 16 | Plotting [10.2s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 28 | saemix parent models [180.2s] +✔ | 28 | saemix parent models [181.5s] ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s] -✔ | 7 | Fitting the SFORB model [3.6s] +✔ | 11 | Processing of residue series +✔ | 7 | Fitting the SFORB model [3.7s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 5 | Summary [0.2s] ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s] -✔ | 9 | Hypothesis tests [7.8s] +✔ | 9 | Hypothesis tests [8.3s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 277.4 s +Duration: 280.7 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1) • Fitting with saemix takes around 10 minutes when using deSolve (1) -[ FAIL 0 | WARN 0 | SKIP 2 | PASS 235 ] +[ FAIL 2 | WARN 0 | SKIP 2 | PASS 244 ] diff --git a/man/set_nd_nq.Rd b/man/set_nd_nq.Rd new file mode 100644 index 00000000..87a3fae1 --- /dev/null +++ b/man/set_nd_nq.Rd @@ -0,0 +1,103 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/set_nd_nq.R +\name{set_nd_nq} +\alias{set_nd_nq} +\alias{set_nd_nq_focus} +\title{Set non-detects and unquantified values in residue series without replicates} +\usage{ +set_nd_nq(res_raw, lod, loq = NA, time_zero_presence = FALSE) + +set_nd_nq_focus( + res_raw, + lod, + loq = NA, + set_first_sample_nd = TRUE, + first_sample_nd_value = 0, + ignore_below_loq_after_first_nd = TRUE +) +} +\arguments{ +\item{res_raw}{Character vector of a residue time series, or matrix of +residue values with rows representing depth profiles for a specific sampling +time, and columns representing time series of residues at the same depth. +Values below the limit of detection (lod) have to be coded as "nd", values +between the limit of detection and the limit of quantification, if any, have +to be coded as "nq". Samples not analysed have to be coded as "na". All +values that are not "na", "nd" or "nq" have to be coercible to numeric} + +\item{lod}{Limit of detection (numeric)} + +\item{loq}{Limit of quantification(numeric). Must be specified if the FOCUS rule to +stop after the first non-detection is to be applied} + +\item{time_zero_presence}{Do we assume that residues occur at time zero? +This only affects samples from the first sampling time that have been +reported as "nd" (not detected).} + +\item{set_first_sample_nd}{Should the first sample be set to "first_sample_nd_value" +in case it is a non-detection?} + +\item{first_sample_nd_value}{Value to be used for the first sample if it is a non-detection} + +\item{ignore_below_loq_after_first_nd}{Should we ignore values below the LOQ after the first +non-detection that occurs after the quantified values?} +} +\value{ +A numeric vector, if a vector was supplied, or a numeric matrix otherwise +} +\description{ +This function automates replacing unquantified values in residue time and +depth series. For time series, the function performs part of the residue +processing proposed in the FOCUS kinetics guidance for parent compounds +and metabolites. For two-dimensional residue series over time and depth, +it automates the proposal of Boesten et al (2015). +} +\section{Functions}{ +\itemize{ +\item \code{set_nd_nq_focus()}: Set non-detects in residue time series according to FOCUS rules + +}} +\examples{ +# FOCUS (2014) p. 75/76 and 131/132 +parent_1 <- c(.12, .09, .05, .03, "nd", "nd", "nd", "nd", "nd", "nd") +set_nd_nq(parent_1, 0.02) +parent_2 <- c(.12, .09, .05, .03, "nd", "nd", .03, "nd", "nd", "nd") +set_nd_nq(parent_2, 0.02) +set_nd_nq_focus(parent_2, 0.02, loq = 0.05) +parent_3 <- c(.12, .09, .05, .03, "nd", "nd", .06, "nd", "nd", "nd") +set_nd_nq(parent_3, 0.02) +set_nd_nq_focus(parent_3, 0.02, loq = 0.05) +metabolite <- c("nd", "nd", "nd", 0.03, 0.06, 0.10, 0.11, 0.10, 0.09, 0.05, 0.03, "nd", "nd") +set_nd_nq(metabolite, 0.02) +set_nd_nq_focus(metabolite, 0.02, 0.05) +# +# Boesten et al. (2015), p. 57/58 +table_8 <- matrix( + c(10, 10, rep("nd", 4), + 10, 10, rep("nq", 2), rep("nd", 2), + 10, 10, 10, "nq", "nd", "nd", + "nq", 10, "nq", rep("nd", 3), + "nd", "nq", "nq", rep("nd", 3), + rep("nd", 6), rep("nd", 6)), + ncol = 6, byrow = TRUE) +set_nd_nq(table_8, 0.5, 1.5, time_zero_presence = TRUE) +table_10 <- matrix( + c(10, 10, rep("nd", 4), + 10, 10, rep("nd", 4), + 10, 10, 10, rep("nd", 3), + "nd", 10, rep("nd", 4), + rep("nd", 18)), + ncol = 6, byrow = TRUE) +set_nd_nq(table_10, 0.5, time_zero_presence = TRUE) +} +\references{ +Boesten, J. J. T. I., van der Linden, A. M. A., Beltman, W. H. +J. and Pol, J. W. (2015). Leaching of plant protection products and their +transformation products; Proposals for improving the assessment of leaching +to groundwater in the Netherlands — Version 2. Alterra report 2630, Alterra +Wageningen UR (University & Research centre) + +FOCUS (2014) Generic Guidance for Estimating Persistence and Degradation +Kinetics from Environmental Fate Studies on Pesticides in EU Registration, Version 1.1, +18 December 2014, p. 251 +} diff --git a/tests/testthat/print_sfo_saem_1.txt b/tests/testthat/print_sfo_saem_1.txt index a11e1e96..34eea058 100644 --- a/tests/testthat/print_sfo_saem_1.txt +++ b/tests/testthat/print_sfo_saem_1.txt @@ -7,13 +7,13 @@ Data: Likelihood computed by importance sampling AIC BIC logLik - 1311 1315 -649 + 1312 1316 -650 Fitted parameters: estimate lower upper -parent_0 1e+02 99.13 1e+02 +parent_0 1e+02 99.28 1e+02 k_parent 4e-02 0.03 4e-02 a.1 9e-01 0.75 1e+00 b.1 5e-02 0.04 5e-02 -SD.parent_0 7e-01 -1.09 3e+00 +SD.parent_0 2e-01 -6.33 7e+00 SD.k_parent 3e-01 0.20 4e-01 diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt index ad701cee..e3e2f7e4 100644 --- a/tests/testthat/summary_hfit_sfo_tc.txt +++ b/tests/testthat/summary_hfit_sfo_tc.txt @@ -45,7 +45,7 @@ SD.log_k_parent 0.3 0.1 0.4 Variance model: est. lower upper -a.1 0.91 0.64 1.17 +a.1 0.90 0.64 1.17 b.1 0.05 0.04 0.06 Backtransformed parameters: diff --git a/tests/testthat/test_set_nd.R b/tests/testthat/test_set_nd.R new file mode 100644 index 00000000..232e7c62 --- /dev/null +++ b/tests/testthat/test_set_nd.R @@ -0,0 +1,86 @@ +context("Processing of residue series") + +# FOCUS (2014) page 76 (parent) and page 132 (metabolite) + +parent_1 <- c(.12, .09, .05, .03, "nd", "nd", "nd", "nd", "nd", "nd") +parent_2 <- c(.12, .09, .05, .03, "nd", "nd", .03, "nd", "nd", "nd") +parent_3 <- c(.12, .09, .05, .03, "nd", "nd", .06, "nd", "nd", "nd") +metabolite <- c("nd", "nd", "nd", 0.03, 0.06, 0.10, 0.11, 0.10, 0.09, 0.05, 0.03, "nd", "nd") + +test_that("Simple residue series are processed as intended", { + + expect_equal(set_nd_nq(parent_1, 0.02), + c(.12, .09, .05, .03, .01, rep(NA, 5))) + + expect_equal(set_nd_nq(parent_2, 0.02, loq = 0.05), + c(.12, .09, .05, .03, .01, .01, .03, .01, NA, NA)) + + expect_equal(set_nd_nq(metabolite, 0.02, loq = 0.05), + c(NA, NA, .01, .03, .06, .1, .11, .1, .09, .05, .03, .01, NA)) + + expect_equal(set_nd_nq(c("nd", 1, 0.2, "nd"), 0.1), + c(NA, 1, 0.2, 0.05)) + +}) + +test_that("Simple residue series are processed as in the FOCUS guidance", { + + # Parent 1 + expect_error(set_nd_nq_focus(parent_1, 0.02), + "You need to specify an LOQ") + expect_equal(set_nd_nq_focus(parent_1, 0.02, 0.05), + c(.12, .09, .05, .03, .01, rep(NA, 5))) + + # Parent 2 + expect_equal(set_nd_nq_focus(parent_2, 0.02, loq = 0.05), + c(.12, .09, .05, .03, .01, rep(NA, 5))) + + # Parent 3 + expect_equal(set_nd_nq_focus(parent_3, 0.02, loq = 0.05), + c(.12, .09, .05, .03, .01, .01, .06, .01, NA, NA)) + + # Metabolite + expect_equal(set_nd_nq_focus(metabolite, 0.02, loq = 0.05), + c(0, NA, .01, .03, .06, .1, .11, .1, .09, .05, .03, .01, NA)) + +}) + +test_that("Examples Boesten et al. (2015, p. 57/58) are correctly processed", { + table_8 <- matrix( + c(10, 10, rep("nd", 4), + 10, 10, rep("nq", 2), rep("nd", 2), + 10, 10, 10, "nq", "nd", "nd", + "nq", 10, "nq", rep("nd", 3), + "nd", "nq", "nq", rep("nd", 3), + rep("nd", 6), rep("nd", 6)), + ncol = 6, byrow = TRUE) + table_8_processed <- set_nd_nq(table_8, 0.5, 1.5, time_zero_presence = TRUE) + table_9 <- matrix( + c(10, 10, 0.25, 0.25, NA, NA, + 10, 10, 1, 1, 0.25, NA, + 10, 10, 10, 1, 0.25, NA, + 1, 10, 1, 0.25, NA, NA, + 0.25, 1, 1, 0.25, NA, NA, + NA, 0.25, 0.25, NA, NA, NA, + rep(NA, 6)), + ncol = 6, byrow = TRUE) + expect_equal(table_8_processed, table_9) + + table_10 <- matrix( + c(10, 10, rep("nd", 4), + 10, 10, rep("nd", 4), + 10, 10, 10, rep("nd", 3), + "nd", 10, rep("nd", 4), + rep("nd", 18)), + ncol = 6, byrow = TRUE) + table_10_processed <- set_nd_nq(table_10, 0.5, time_zero_presence = TRUE) + table_11 <- matrix( + c(10, 10, 0.25, rep(NA, 3), + 10, 10, 0.25, rep(NA, 3), + 10, 10, 10, 0.25, NA, NA, + 0.25, 10, 0.25, rep(NA, 3), + NA, 0.25, rep(NA, 4), + rep(NA, 12)), + ncol = 6, byrow = TRUE) + expect_equal(table_10_processed, table_11) +}) |