From b847fec686bc1db59079412eb18063d3514ecf75 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 21 Dec 2016 10:35:16 +0100 Subject: TWA calculation for SFO, FOMC and DFOP (parent) --- DESCRIPTION | 4 +-- NEWS.md | 6 +++++ R/twa.R | 65 +++++++++++++++++++++++++++++++++++++++++++++++ build.log | 7 +---- man/twa.Rd | 40 +++++++++++++++++++++++++++++ test.log | 1 + tests/testthat/test_twa.R | 45 ++++++++++++++++++++++++++++++++ 7 files changed, 160 insertions(+), 8 deletions(-) create mode 100644 R/twa.R create mode 100644 man/twa.Rd create mode 100644 tests/testthat/test_twa.R diff --git a/DESCRIPTION b/DESCRIPTION index 55467688..cdc1a4a5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data -Version: 0.9.45 -Date: 2016-12-08 +Version: 0.9.45.1 +Date: 2016-12-21 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de"), person("Katrin", "Lindenberger", role = "ctb"), diff --git a/NEWS.md b/NEWS.md index 5397dfce..ac36f39b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# mkin 0.9.45.1 (2016-12-20) + +## New features + +- A `twa` function, calculating maximum time weighted average concentrations for the parent (SFO, FOMC and DFOP). + # mkin 0.9.45 (2016-12-08) ## Minor changes diff --git a/R/twa.R b/R/twa.R new file mode 100644 index 00000000..e8894b0b --- /dev/null +++ b/R/twa.R @@ -0,0 +1,65 @@ +# Copyright (C) 2016 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 + +twa <- function(fit, windows) { + parms.all <- c(fit$bparms.optim, fit$bparms.fixed) + obs_vars <- fit$obs_vars + if (length(obs_vars) > 1) { + warning("Calculation of maximum time weighted average concentrations is", + "currently only implemented for the parent compound using", + "analytical solutions") + } + obs_var <- obs_vars[1] + spec = fit$mkinmod$spec + type = spec[[1]]$type + + M0 <- parms.all[paste0(obs_var, "_0")] + + if (type == "SFO") { + k_name <- paste0("k_", obs_var) + if (fit$mkinmod$use_of_ff == "min") { + k_name <- paste0(k_name, "_sink") + } + k <- parms.all[k_name] + twafunc <- function(t) { + M0 * (1 - exp(- k * t)) / (k * t) + } + } + if (type == "FOMC") { + alpha <- parms.all["alpha"] + beta <- parms.all["beta"] + twafunc <- function(t) { + M0 * (beta)/(t * (1 - alpha)) * ((t/beta + 1)^(1 - alpha) - 1) + } + } + if (type == "DFOP") { + k1 <- parms.all["k1"] + k2 <- parms.all["k2"] + g <- parms.all["g"] + twafunc <- function(t) { + M0/t * ((g/k1) * (1 - exp(- k1 * t)) + ((1 - g)/k2) * (1 - exp(- k2 * t))) + } + } + if (type %in% c("HS", "IORE", "SFORB")) { + stop("Calculation of maximum time weighted average concentrations is currently ", + "not implemented for the ", type, " model.") + } + res <- twafunc(windows) + names(res) <- windows + return(res) +} diff --git a/build.log b/build.log index efdaa3b6..a02e93ae 100644 --- a/build.log +++ b/build.log @@ -2,9 +2,4 @@ * preparing ‘mkin’: * checking DESCRIPTION meta-information ... OK * installing the package to build vignettes -* creating vignettes ... OK -* checking for LF line-endings in source and make files -* checking for empty or unneeded directories -* looking to see if a ‘data/datalist’ file should be added -* building ‘mkin_0.9.45.tar.gz’ - +* creating vignettes ... \ No newline at end of file diff --git a/man/twa.Rd b/man/twa.Rd new file mode 100644 index 00000000..57413b18 --- /dev/null +++ b/man/twa.Rd @@ -0,0 +1,40 @@ +\name{twa} +\alias{twa} +\title{ + Function to calculate maximum time weighted average concentrations from kinetic models fitted with mkinfit +} +\description{ +This function calculates maximum moving window time weighted average concentrations +(TWAs) for kinetic models fitted with \code{\link{mkinfit}}. Currently, only +calculations for the parent are implemented for the SFO, FOMC and DFOP models, +using the analytical formulas given in the PEC soil section of the FOCUS +guidance.} +\usage{ + twa(fit, windows) +} +\arguments{ + \item{fit}{ + An object of class \code{\link{mkinfit}}. + } + \item{windows}{ + The width of the time windows for which the TWAs should be calculated. + } +} +\value{ + A numeric vector, named using the \code{windows} argument. +} +\examples{ + fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) + twa(fit, c(7, 21)) +} +\references{ + FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and + Degradation Kinetics from Environmental Fate Studies on Pesticides in EU + Registration} Report of the FOCUS Work Group on Degradation Kinetics, + EC Document Reference Sanco/10058/2005 version 2.0, 434 pp, + \url{http://http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics} +} +\author{ + Johannes Ranke +} +\keyword{ manip } diff --git a/test.log b/test.log index 436b9663..935df857 100644 --- a/test.log +++ b/test.log @@ -14,6 +14,7 @@ Model predictions with mkinpredict: ... Fitting of parent only models: ..................... Complex test case from Schaefer et al. (2007) Piacenza paper: .. Results for synthetic data established in expertise for UBA (Ranke 2014): .... +Calculation of maximum time weighted average concentrations (TWAs): ... Skipped ------------------------------------------------------------------------ 1. Fitting with large parameter correlation gives warnings (@test_FOMC_ill-defined.R#30) - Skip test for warnings triggered by large parameter correlation as it failed on r-forge diff --git a/tests/testthat/test_twa.R b/tests/testthat/test_twa.R new file mode 100644 index 00000000..3a28e0c1 --- /dev/null +++ b/tests/testthat/test_twa.R @@ -0,0 +1,45 @@ +# Copyright (C) 2016 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 + +context("Calculation of maximum time weighted average concentrations (TWAs)") + +twa_models <- c("SFO", "FOMC", "DFOP") +fits <- mmkin(twa_models, list(FOCUS_D = FOCUS_2006_D), + quiet = TRUE) + +SFO_SFO <- mkinmod(parent = list(type = "SFO", to = "m1"), + m1 = list(type = "SFO"), quiet = TRUE) +fit.m1 <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) + +test_that("Time weighted average concentrations are correct", { + outtimes_7 <- seq(0, 7, length.out = 10000) + for (model in twa_models) { + fit <- fits[[model, 1]] + bpar <- summary(fit)$bpar[, "Estimate"] + pred_7 <- mkinpredict(fit$mkinmod, + odeparms = bpar[2:length(bpar)], + odeini = c(parent = bpar[[1]]), + outtimes = outtimes_7) + twa_num <- mean(pred_7$parent) + names(twa_num) <- 7 + twa_ana <- twa(fit, 7) + + # Test for absolute difference (scale = 1) + expect_equal(twa_num, twa_ana, tolerance = 0.001, scale = 1) + } +}) -- cgit v1.2.1