diff options
| -rw-r--r-- | .Rbuildignore | 1 | ||||
| -rw-r--r-- | ChangeLog | 6 | ||||
| -rw-r--r-- | NAMESPACE | 3 | ||||
| -rw-r--r-- | R/twa.R | 88 | ||||
| -rw-r--r-- | docs/reference/max_twa.html | 45 | ||||
| -rw-r--r-- | docs/reference/one_box-4.png | bin | 5743 -> 5587 bytes | |||
| -rw-r--r-- | docs/reference/one_box-8.png | bin | 9901 -> 9860 bytes | |||
| -rw-r--r-- | docs/reference/one_box.html | 12 | ||||
| -rw-r--r-- | docs/reference/plot.one_box-8.png | bin | 13603 -> 13490 bytes | |||
| -rw-r--r-- | docs/reference/sawtooth-6.png | bin | 13603 -> 13490 bytes | |||
| -rw-r--r-- | docs/reference/sawtooth.html | 4 | ||||
| -rw-r--r-- | man/max_twa.Rd | 19 | ||||
| -rw-r--r-- | man/one_box.Rd | 12 | ||||
| -rw-r--r-- | tests/testthat/test_max_twa.R | 11 |
14 files changed, 182 insertions, 19 deletions
diff --git a/.Rbuildignore b/.Rbuildignore index 48acd64..59d7984 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,4 +7,5 @@ ^build.log$ ^test.log$ ^inst/extdata/Tabelle\ der\ Abdrifteckwerte.xls$ +^_pkgdown.yml$ ^docs$ @@ -1,3 +1,9 @@ +commit 3ead7acba845b4f2552f555dfb29da889ed0cda8 +Author: Johannes Ranke <jranke@uni-bremen.de> +Date: 2017-01-19 09:42:21 +0100 + + Make max_twa() a bit safer + commit 74ed85b07f09ea99476208749cd274a476ba4536 Author: Johannes Ranke <jranke@uni-bremen.de> Date: 2017-01-19 09:27:36 +0100 @@ -2,6 +2,7 @@ S3method(GUS,chent) S3method(GUS,numeric) +S3method(max_twa,mkinfit) S3method(max_twa,one_box) S3method(one_box,character) S3method(one_box,mkinfit) @@ -37,8 +38,10 @@ importFrom(methods,is) importFrom(mkin,mkinpredict) importFrom(readr,fwf_empty) importFrom(readr,read_fwf) +importFrom(stats,end) importFrom(stats,filter) importFrom(stats,frequency) importFrom(stats,plot.ts) +importFrom(stats,start) importFrom(stats,time) importFrom(stats,ts) @@ -1,4 +1,4 @@ -# Copyright (C) 2017 Johannes Ranke +# Copyright (C) 2016,2017 Johannes Ranke # Contact: jranke@uni-bremen.de # This file is part of the R package pfm @@ -23,6 +23,7 @@ #' #' @param x When numeric, this is the half-life to be used for an exponential #' decline. If x is an mkinfit object, the decline is calculated from this object +#' @param ini The initial amount for each compound #' @param t_end End of the time series #' @param res Resolution of the time series #' @param ... Further arguments passed to methods @@ -44,7 +45,7 @@ #' fit_2 <- mkinfit(m_2, FOCUS_2006_D, quiet = TRUE) #' pred_2 <- one_box(fit_2) #' plot(pred_2) -one_box <- function(x, ..., +one_box <- function(x, ini, ..., t_end = 100, res = 0.01) { UseMethod("one_box") @@ -52,7 +53,7 @@ one_box <- function(x, ..., #' @rdname one_box #' @export -one_box.numeric <- function(x, ..., +one_box.numeric <- function(x, ini = 1, ..., t_end = 100, res = 0.01) { half_life = x @@ -68,7 +69,7 @@ one_box.numeric <- function(x, ..., #' @rdname one_box #' @param parms A named numeric vector containing the model parameters #' @export -one_box.character <- function(x, parms, ..., +one_box.character <- function(x, ini = 1, parms, ..., t_end = 100, res = 0.01) { parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE") @@ -96,11 +97,18 @@ one_box.character <- function(x, parms, ..., #' @rdname one_box #' @importFrom mkin mkinpredict #' @export -one_box.mkinfit <- function(x, ..., t_end = 100, res = 0.01) { +one_box.mkinfit <- function(x, ini = c("model", 1), ..., t_end = 100, res = 0.01) { fit <- x + + ini = match.arg(ini) + if (ini == "model") { + odeini = x$bparms.state + } else { + odeini <- c(1, rep(0, length(fit$mkinmod$spec) - 1)) + names(odeini) <- names(fit$mkinmod$spec) + } + t_out = seq(0, t_end, by = res) - odeini <- c(1, rep(0, length(fit$mkinmod$spec) - 1)) - names(odeini) <- names(fit$mkinmod$spec) if (length(fit$mkinmod$spec) == 1) solution_type = "analytical" else solution_type = "deSolve" @@ -216,6 +224,7 @@ sawtooth <- function(x, n = 1, i = 365, #' @param x An object of type \code{\link{one_box}} #' @param window The size of the moving window #' @seealso \code{\link{max_twa}} +#' @importFrom stats start end #' @export #' @examples #' pred <- sawtooth(one_box(10), @@ -247,12 +256,77 @@ twa.one_box <- function(x, window = 21) #' for finding the maximum. It is therefore recommended to check this using #' \code{\link{plot.one_box}} using the window size for the argument #' \code{max_twa}. +#' +#' The method working directly on fitted \code{\link{mkinfit}} objects uses the +#' equations given in the PEC soil section of the FOCUS guidance and is restricted +#' SFO, FOMC and DFOP models and to the parent compound +#' @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://esdac.jrc.ec.europa.eu/projects/degradation-kinetics} #' @seealso \code{\link{twa}} #' @inheritParams twa #' @export +#' @examples +#' pred <- sawtooth(one_box(10), +#' applications = data.frame(time = c(0, 7), amount = c(1, 1))) +#' max_twa(pred) +#' pred_FOMC <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) +#' max_twa(pred_FOMC) max_twa <- function(x, window = 21) UseMethod("max_twa") #' @export +max_twa.mkinfit <- function(x, window = 21) { + fit <- x + 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(t = window) + names(res) <- window + return(res) +} + +#' @export max_twa.one_box <- function(x, window = 21) { freq = frequency(x) diff --git a/docs/reference/max_twa.html b/docs/reference/max_twa.html index 4bb1be9..f78da58 100644 --- a/docs/reference/max_twa.html +++ b/docs/reference/max_twa.html @@ -70,7 +70,11 @@ </div> - <p>The maximum time weighted average concentration for a moving window</p> + <p>If you generate your time series using <code><a href='sawtooth.html'>sawtooth</a></code>, +you need to make sure that the length of the time series allows +for finding the maximum. It is therefore recommended to check this using +<code><a href='plot.one_box.html'>plot.one_box</a></code> using the window size for the argument +<code>max_twa</code>.</p> <pre><span class='fu'>max_twa</span>(<span class='no'>x</span>, <span class='kw'>window</span> <span class='kw'>=</span> <span class='fl'>21</span>)</pre> @@ -83,19 +87,56 @@ <dd>The size of the moving window</dd> </dl> + <h2 class="hasAnchor" id="details"><a class="anchor" href="#details"></a>Details</h2> + + <p>The method working directly on fitted <code>mkinfit</code> objects uses the +equations given in the PEC soil section of the FOCUS guidance and is restricted +SFO, FOMC and DFOP models and to the parent compound</p> + + <h2 class="hasAnchor" id="references"><a class="anchor" href="#references"></a>References</h2> + + <p>FOCUS (2006) “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, + <a href = 'http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics'>http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics</a></p> + <h2 class="hasAnchor" id="see-also"><a class="anchor" href="#see-also"></a>See also</h2> <p><code><a href='twa.html'>twa</a></code></p> + <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> + <pre class="examples"><div class='input'><span class='no'>pred</span> <span class='kw'><-</span> <span class='fu'><a href='sawtooth.html'>sawtooth</a></span>(<span class='fu'><a href='one_box.html'>one_box</a></span>(<span class='fl'>10</span>), + <span class='kw'>applications</span> <span class='kw'>=</span> <span class='fu'>data.frame</span>(<span class='kw'>time</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='fl'>0</span>, <span class='fl'>7</span>), <span class='kw'>amount</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='fl'>1</span>, <span class='fl'>1</span>))) +<span class='fu'>max_twa</span>(<span class='no'>pred</span>)</div><div class='output co'>#> $max +#> parent +#> 0.9537545 +#> +#> $window_start +#> parent +#> 0 +#> +#> $window_end +#> parent +#> 21 +#> </div><div class='input'><span class='no'>pred_FOMC</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='st'>"FOMC"</span>, <span class='no'>FOCUS_2006_C</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) +<span class='fu'>max_twa</span>(<span class='no'>pred_FOMC</span>)</div><div class='output co'>#> 21 +#> 18.22124 </div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> <h2>Contents</h2> <ul class="nav nav-pills nav-stacked"> <li><a href="#arguments">Arguments</a></li> + <li><a href="#details">Details</a></li> + + <li><a href="#references">References</a></li> + <li><a href="#see-also">See also</a></li> - </ul> + + <li><a href="#examples">Examples</a></li> + </ul> </div> </div> diff --git a/docs/reference/one_box-4.png b/docs/reference/one_box-4.png Binary files differindex 10ebb35..3c6dadb 100644 --- a/docs/reference/one_box-4.png +++ b/docs/reference/one_box-4.png diff --git a/docs/reference/one_box-8.png b/docs/reference/one_box-8.png Binary files differindex 31e2f4b..897737e 100644 --- a/docs/reference/one_box-8.png +++ b/docs/reference/one_box-8.png diff --git a/docs/reference/one_box.html b/docs/reference/one_box.html index 698f887..cb49c02 100644 --- a/docs/reference/one_box.html +++ b/docs/reference/one_box.html @@ -74,22 +74,26 @@ This does not create objects of type <code>ts</code>.</p> - <pre><span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) + <pre><span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='no'>ini</span>, <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) <span class='co'># S3 method for numeric</span> -<span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) +<span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='kw'>ini</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) <span class='co'># S3 method for character</span> -<span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='no'>parms</span>, <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) +<span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='kw'>ini</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='no'>parms</span>, <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, + <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) <span class='co'># S3 method for mkinfit</span> -<span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>)</pre> +<span class='fu'>one_box</span>(<span class='no'>x</span>, <span class='kw'>ini</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='st'>"model"</span>, <span class='fl'>1</span>), <span class='no'>...</span>, <span class='kw'>t_end</span> <span class='kw'>=</span> <span class='fl'>100</span>, + <span class='kw'>res</span> <span class='kw'>=</span> <span class='fl'>0.01</span>)</pre> <h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a> Arguments</h2> <dl class="dl-horizontal"> <dt>x</dt> <dd>When numeric, this is the half-life to be used for an exponential decline. If x is an mkinfit object, the decline is calculated from this object</dd> + <dt>ini</dt> + <dd>The initial amount for each compound</dd> <dt>...</dt> <dd>Further arguments passed to methods</dd> <dt>t_end</dt> diff --git a/docs/reference/plot.one_box-8.png b/docs/reference/plot.one_box-8.png Binary files differindex 161ecb4..f534444 100644 --- a/docs/reference/plot.one_box-8.png +++ b/docs/reference/plot.one_box-8.png diff --git a/docs/reference/sawtooth-6.png b/docs/reference/sawtooth-6.png Binary files differindex 161ecb4..f534444 100644 --- a/docs/reference/sawtooth-6.png +++ b/docs/reference/sawtooth-6.png diff --git a/docs/reference/sawtooth.html b/docs/reference/sawtooth.html index cf33630..216f64a 100644 --- a/docs/reference/sawtooth.html +++ b/docs/reference/sawtooth.html @@ -102,8 +102,8 @@ the corresponding amounts applied in the second column.</dd> <span class='no'>pred_2_saw</span> <span class='kw'><-</span> <span class='fu'>sawtooth</span>(<span class='no'>pred_2</span>, <span class='fl'>2</span>, <span class='fl'>7</span>) <span class='fu'>plot</span>(<span class='no'>pred_2_saw</span>, <span class='kw'>max_twa</span> <span class='kw'>=</span> <span class='fl'>21</span>, <span class='kw'>max_twa_var</span> <span class='kw'>=</span> <span class='st'>"m1"</span>)</div><img src='sawtooth-6.png' alt='' width='540' height='400' /><div class='input'> <span class='fu'><a href='max_twa.html'>max_twa</a></span>(<span class='no'>pred_2_saw</span>)</div><div class='output co'>#> $max -#> parent m1 -#> 0.6627707 0.8542785 +#> parent m1 +#> 66.01096 85.08484 #> #> $window_start #> parent m1 diff --git a/man/max_twa.Rd b/man/max_twa.Rd index e39165a..0a8dcca 100644 --- a/man/max_twa.Rd +++ b/man/max_twa.Rd @@ -18,6 +18,25 @@ for finding the maximum. It is therefore recommended to check this using \code{\link{plot.one_box}} using the window size for the argument \code{max_twa}. } +\details{ +The method working directly on fitted \code{\link{mkinfit}} objects uses the +equations given in the PEC soil section of the FOCUS guidance and is restricted +SFO, FOMC and DFOP models and to the parent compound +} +\examples{ +pred <- sawtooth(one_box(10), + applications = data.frame(time = c(0, 7), amount = c(1, 1))) +max_twa(pred) +pred_FOMC <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) +max_twa(pred_FOMC) +} +\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://esdac.jrc.ec.europa.eu/projects/degradation-kinetics} +} \seealso{ \code{\link{twa}} } diff --git a/man/one_box.Rd b/man/one_box.Rd index 5a46142..8132a70 100644 --- a/man/one_box.Rd +++ b/man/one_box.Rd @@ -7,18 +7,22 @@ \alias{one_box.mkinfit} \title{Create a time series of decline data} \usage{ -one_box(x, ..., t_end = 100, res = 0.01) +one_box(x, ini, ..., t_end = 100, res = 0.01) -\method{one_box}{numeric}(x, ..., t_end = 100, res = 0.01) +\method{one_box}{numeric}(x, ini = 1, ..., t_end = 100, res = 0.01) -\method{one_box}{character}(x, parms, ..., t_end = 100, res = 0.01) +\method{one_box}{character}(x, ini = 1, parms, ..., t_end = 100, + res = 0.01) -\method{one_box}{mkinfit}(x, ..., t_end = 100, res = 0.01) +\method{one_box}{mkinfit}(x, ini = c("model", 1), ..., t_end = 100, + res = 0.01) } \arguments{ \item{x}{When numeric, this is the half-life to be used for an exponential decline. If x is an mkinfit object, the decline is calculated from this object} +\item{ini}{The initial amount for each compound} + \item{...}{Further arguments passed to methods} \item{t_end}{End of the time series} diff --git a/tests/testthat/test_max_twa.R b/tests/testthat/test_max_twa.R new file mode 100644 index 0000000..82d20f7 --- /dev/null +++ b/tests/testthat/test_max_twa.R @@ -0,0 +1,11 @@ +library(pfm) +context("Check max_twa for parent mkinfit models against analytical solutions") + +test_that("max_twa for simple decline curves is correct", { + fits <- mmkin(c("SFO", "FOMC", "DFOP"), list(FOCUS_2006_C), cores = 1, quiet = TRUE) + max_twa_21_analytical <- sapply(fits, max_twa, 21) + names(max_twa_21_analytical) <- rep("parent", 3) + preds <- lapply(fits, one_box) + max_twa_21_preds <- sapply(preds, function(x) max_twa(x, 21)[["max"]]) + expect_equal(max_twa_21_analytical, max_twa_21_preds, tolerance = .03, scale = 1) +}) |
