aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2017-01-19 10:36:19 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2017-01-19 10:36:19 +0100
commitdd30f0d0ff1d8d0cc46aaef6e0917c51fe798f52 (patch)
treea954197c4e8e093d4b08c08eab14753b58602ff4
parent3ead7acba845b4f2552f555dfb29da889ed0cda8 (diff)
Move mkin::twa to pfm::max_twa.mkinfit
- Add max_twa.mkinfit() recently introduced to mkin as mkin::twa() but never released with it - Add a test to check max_twa.one_box() against analytical solutions in max_twa.mkinfit(). - Clean up R CMD check - Update docs
-rw-r--r--.Rbuildignore1
-rw-r--r--ChangeLog6
-rw-r--r--NAMESPACE3
-rw-r--r--R/twa.R88
-rw-r--r--docs/reference/max_twa.html45
-rw-r--r--docs/reference/one_box-4.pngbin5743 -> 5587 bytes
-rw-r--r--docs/reference/one_box-8.pngbin9901 -> 9860 bytes
-rw-r--r--docs/reference/one_box.html12
-rw-r--r--docs/reference/plot.one_box-8.pngbin13603 -> 13490 bytes
-rw-r--r--docs/reference/sawtooth-6.pngbin13603 -> 13490 bytes
-rw-r--r--docs/reference/sawtooth.html4
-rw-r--r--man/max_twa.Rd19
-rw-r--r--man/one_box.Rd12
-rw-r--r--tests/testthat/test_max_twa.R11
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$
diff --git a/ChangeLog b/ChangeLog
index 3aff70c..f7a68ea 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -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
diff --git a/NAMESPACE b/NAMESPACE
index c4d671f..35245d9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -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)
diff --git a/R/twa.R b/R/twa.R
index fcf49c0..88cfec3 100644
--- a/R/twa.R
+++ b/R/twa.R
@@ -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) &#8220;Guidance Document on Estimating Persistence and
+ Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
+ Registration&#8221; 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'>&lt;-</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'>#&gt; $max
+#&gt; parent
+#&gt; 0.9537545
+#&gt;
+#&gt; $window_start
+#&gt; parent
+#&gt; 0
+#&gt;
+#&gt; $window_end
+#&gt; parent
+#&gt; 21
+#&gt; </div><div class='input'><span class='no'>pred_FOMC</span> <span class='kw'>&lt;-</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'>#&gt; 21
+#&gt; 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
index 10ebb35..3c6dadb 100644
--- a/docs/reference/one_box-4.png
+++ b/docs/reference/one_box-4.png
Binary files differ
diff --git a/docs/reference/one_box-8.png b/docs/reference/one_box-8.png
index 31e2f4b..897737e 100644
--- a/docs/reference/one_box-8.png
+++ b/docs/reference/one_box-8.png
Binary files differ
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
index 161ecb4..f534444 100644
--- a/docs/reference/plot.one_box-8.png
+++ b/docs/reference/plot.one_box-8.png
Binary files differ
diff --git a/docs/reference/sawtooth-6.png b/docs/reference/sawtooth-6.png
index 161ecb4..f534444 100644
--- a/docs/reference/sawtooth-6.png
+++ b/docs/reference/sawtooth-6.png
Binary files differ
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'>&lt;-</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'>#&gt; $max
-#&gt; parent m1
-#&gt; 0.6627707 0.8542785
+#&gt; parent m1
+#&gt; 66.01096 85.08484
#&gt;
#&gt; $window_start
#&gt; 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)
+})

Contact - Imprint