TWA calculation for SFO, FOMC and DFOP (parent)

7 files changed, 160 insertions, 8 deletions

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 <http://www.gnu.org/licenses/>
+
+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 <http://www.gnu.org/licenses/>
+
+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)
+  }
+})