mkindsg class to hold groups of datasets

- D24_2014 dataset on aerobic soil degradation of 2,4-D from the EU
  assessment as mkindsg object with metadata
- f_time_norm_focus() to do time-step normalisation using the FOCUS
  method
- focus_soil_moisture data with default moisture contents at pF1,
  pF 2 and pF 2.5 for USDA soil types from FOCUS GW guidance
- Dataset generation scripts in inst/dataset_generation
- Depend on R >= 2.15.1 in order to facilitate the use of
  utils::globalVariables()

12 files changed, 233 insertions, 16 deletions

diff --git a/R/D24_2014.R b/R/D24_2014.R
new file mode 100644
index 00000000..81316833
--- /dev/null
+++ b/R/D24_2014.R
@@ -0,0 +1,28 @@
+#'  Aerobic soil degradation data on 2,4-D from the EU assessment in 2014
+#'
+#' The five datasets were extracted from the active substance evaluation dossier
+#' published by EFSA. Kinetic evaluations shown for these datasets are intended
+#' to illustrate and advance kinetic modelling. The fact that these data and
+#' some results are shown here does not imply a license to use them in the
+#' context of pesticide  registrations, as the use of the data may be
+#' constrained by data protection regulations.
+#'
+#' Metabolite residues at early sampling times reported as 0.0 were set to NA.
+#'
+#' The R code used to create this data object is installed with this package
+#' in the 'dataset_generation' directory. In the code, page numbers are given for
+#' specific pieces of information in the comments.
+#'
+#' @format An [mkindsg] object grouping five datasets
+#' @source Hellenic Ministry of Rural Development and Agriculture (2014)
+#'   Final addendum to the Renewal Assessment Report - public version - 2,4-D
+#'   Volume 3 Annex B.8 Fate and behaviour in the environment p. 638, 640,
+#'   644-646.
+#'   \url{http://registerofquestions.efsa.europa.eu/roqFrontend/outputLoader?output=ON-3812}
+#' @examples
+#' print(D24_2014)
+#' print(D24_2014$ds[[1]], data = TRUE)
+#' m1 = mkinmod(D24 = list(type = "SFO", to = "phenol"),
+#'   phenol = list(type = "SFO", to = "anisole"),
+#'   anisole = list(type = "SFO"))
+"D24_2014"
diff --git a/R/f_time_norm_focus.R b/R/f_time_norm_focus.R
new file mode 100644
index 00000000..be5cf583
--- /dev/null
+++ b/R/f_time_norm_focus.R
@@ -0,0 +1,77 @@
+utils::globalVariables("D24_2014")
+
+#' Normalisation factors for aerobic soil degradation according to FOCUS guidance
+#'
+#' Time step normalisation factors for aerobic soil degradation as described
+#' in Appendix 8 to the FOCUS kinetics guidance (FOCUS 2014, p. 369).
+#'
+#' @param object An object containing information used for the calculations
+#' @param temperature Numeric vector of temperatures in °C
+#' @param moisture Numeric vector of moisture contents in \\% w/w
+#' @param field_moisture Numeric vector of moisture contents at field capacity
+#'   (pF2) in \\% w/w
+#' @param Q10 The Q10 value used for temperature normalisation
+#' @param walker The Walker exponent used for moisture normalisation
+#' @param f_na The factor to use for NA values. If set to NA, only factors
+#'  for complete cases will be returned.
+#' @param \dots Currently not used
+#' @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}
+#' FOCUS (2014) \dQuote{Generic guidance for Estimating Persistence
+#'   and Degradation Kinetics from Environmental Fate Studies on Pesticides in
+#'   EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
+#'   Version 1.1, 18 December 2014
+#'   \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
+#' @seealso [focus_soil_moisture]
+#' @examples
+#' f_time_norm_focus(25, 20, 25) # 1.37, compare p. 184
+#'
+#' D24_2014$meta
+#' # No moisture normalisation in the first dataset, so we use f_na = 1 to get
+#' # Temperature only normalisation as in the EU evaluation
+#' f_time_norm_focus(D24_2014, f_na = 1)
+#' # Moisture normalisation for the other four soils is one, as soil moisture
+#' # is higher than the approximate field capacity derived from the USDA soil
+#' # type
+#' @export
+f_time_norm_focus <- function(object, ...) {
+  UseMethod("f_time_norm_focus")
+}
+
+#' @rdname f_time_norm_focus
+#' @export
+f_time_norm_focus.numeric <- function(object, 
+  moisture = NA, field_moisture = NA, 
+  temperature = object, 
+  Q10 = 2.58, walker = 0.7, f_na = NA, ...)
+{
+  f_temp <- ifelse(is.na(temperature),
+    f_na,
+    ifelse(temperature <= 0,
+      0,
+      Q10^((temperature - 20)/10)))
+  f_moist <- ifelse(is.na(moisture),
+    f_na,
+    ifelse(moisture >= field_moisture,
+      1,
+      (moisture / field_moisture)^walker))
+  f_time_norm <- f_temp * f_moist
+  f_time_norm
+}
+
+#' @rdname f_time_norm_focus
+#' @export
+f_time_norm_focus.mkindsg <- function(object, Q10 = 2.58, walker = 0.7, f_na = NA, ...) {
+  meta <- object$meta
+  field_moisture <- focus_soil_moisture[meta$usda_soil_type, "pF2"]
+  study_moisture <- meta$rel_moisture * meta$moisture_ref
+  object$f_time_norm <- f_time_norm_focus(meta$temperature,
+    moisture = study_moisture, field_moisture = field_moisture,
+    Q10 = Q10, walker = walker, f_na = f_na)
+  cat("$time_norm was set to\n")
+  print(object$f_time_norm)
+}
diff --git a/R/focus_soil_moisture.R b/R/focus_soil_moisture.R
new file mode 100644
index 00000000..7b22fbcc
--- /dev/null
+++ b/R/focus_soil_moisture.R
@@ -0,0 +1,11 @@
+utils::globalVariables("focus_soil_moisture")
+
+#' FOCUS default values for soil moisture contents at field capacity, MWHC and 1/3 bar
+#'
+#' @format A matrix with upper case USDA soil classes as row names, and water tension
+#'   ('pF1', 'pF2', 'pF 2.5') as column names
+#' @source Anonymous (2014) Generic Guidance for Tier 1 FOCUS Ground Water Assessment
+#'   Version 2.2, May 2014 \url{https://esdac.jrc.ec.europa.eu/projects/ground-water}
+#' @examples
+#' focus_soil_moisture
+"focus_soil_moisture"
diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R
index 971f5273..3d0428d2 100644
--- a/R/mkin_wide_to_long.R
+++ b/R/mkin_wide_to_long.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
+utils::globalVariables(c("name", "time", "value"))
 
 #' Convert a dataframe with observations over time into long format
 #' 
diff --git a/R/mkinds.R b/R/mkinds.R
index d6f296bf..0e970694 100644
--- a/R/mkinds.R
+++ b/R/mkinds.R
@@ -1,5 +1,5 @@
 #' A dataset class for mkin
-#' 
+#'
 #' @description
 #' At the moment this dataset class is hardly used in mkin. For example,
 #' mkinfit does not take mkinds datasets as argument, but works with dataframes
@@ -7,12 +7,11 @@
 #' provided by this package come as mkinds objects nevertheless.
 #'
 #' @importFrom R6 R6Class
-#' @seealso The S3 printing method \code{\link{print.mkinds}}
 #' @examples
-#' 
+#'
 #' mds <- mkinds$new("FOCUS A", FOCUS_2006_A)
 #' print(mds)
-#' 
+#'
 #' @export
 mkinds <- R6Class("mkinds",
   public = list(
@@ -62,15 +61,117 @@ mkinds <- R6Class("mkinds",
 )
 
 #' Print mkinds objects
-#' 
-#' @param x An \code{\link{mkinds}} object.
+#'
+#' @rdname mkinds
+#' @param x An [mkinds] object.
+#' @param data Should the data be printed?
 #' @param \dots Not used.
 #' @export
-print.mkinds <- function(x, ...) {
+print.mkinds <- function(x, data = FALSE, ...) {
   cat("<mkinds> with $title: ",  x$title, "\n")
   cat("Observed compounds $observed: ", paste(x$observed, collapse = ", "), "\n")
-  cat("Sampling times $sampling_times: ", paste(x$sampling_times, collapse = ", "), "\n")
+  cat("Sampling times $sampling_times:\n")
+  cat(paste(x$sampling_times, collapse = ", "), "\n")
   cat("With a maximum of ", x$replicates, " replicates\n")
   if (!is.na(x$time_unit)) cat("Time unit: ", x$time_unit, "\n")
   if (!is.na(x$unit)) cat("Observation unit: ", x$unit, "\n")
+  if (data) print(mkin_long_to_wide(x$data))
+}
+
+#' A class for dataset groups for mkin
+#'
+#' @description
+#' A container for working with datasets that share at least one compound,
+#' so that combined evaluations are desirable.
+#'
+#' Time normalisation factors are initialised with a value of 1 for each
+#' dataset if no data are supplied.
+#'
+#' @examples
+#'
+#' mdsg <- mkindsg$new("Experimental X", experimental_data_for_UBA_2019[6:10])
+#' print(mdsg)
+#' print(mdsg, verbose = TRUE)
+#' print(mdsg, verbose = TRUE, data = TRUE)
+#'
+#' @export
+mkindsg <- R6Class("mkindsg",
+  public = list(
+
+    #' @field title A title for the dataset group
+    title = NULL,
+
+    #' @field ds A list of mkinds objects
+    ds = NULL,
+
+    #' @field observed_n Occurrence counts of compounds in datasets
+    observed_n = NULL,
+
+    #' @field f_time_norm Time normalisation factors
+    f_time_norm = NULL,
+
+    #' @field meta A data frame with a row for each dataset,
+    #'   containing additional information in the form
+    #'   of categorical data (factors) or numerical data
+    #'   (e.g. temperature, moisture,
+    #'   or covariates like soil pH).
+    meta = NULL,
+
+    #' @description
+    #' Create a new mkindsg object
+    #' @param title The title
+    #' @param ds A list of mkinds objects
+    #' @param f_time_norm Time normalisation factors
+    #' @param meta The meta data
+    initialize = function(title = "", ds,
+      f_time_norm = rep(1, length(ds)), meta)
+    {
+      self$title <- title
+      if (all(sapply(ds, inherits, "mkinds"))) {
+        self$ds <- ds
+      } else {
+        stop("Please supply a list of mkinds objects")
+      }
+
+      all_observed <- unlist(lapply(ds, function(x) x$observed))
+      observed <- factor(all_observed, levels = unique(all_observed))
+      self$observed_n <- table(observed)
+      names(dimnames(self$observed_n)) <- NULL
+      self$f_time_norm <- f_time_norm
+
+      if (!missing(meta)) {
+        self$meta <- meta
+      }
+    }
+  )
+)
+
+#' Print mkindsg objects
+#'
+#' @rdname mkindsg
+#' @param x An [mkindsg] object.
+#' @param verbose Should the mkinds objects be printed?
+#' @param data Should the mkinds objects be printed with their data?
+#' @param \dots Not used.
+#' @export
+print.mkindsg <- function(x, data = FALSE, verbose = data, ...) {
+  cat("<mkindsg> holding", length(x$ds), "mkinds objects\n")
+  cat("Title $title: ",  x$title, "\n")
+  cat("Occurrene of observed compounds $observed_n:\n")
+  print(x$observed_n)
+  if (any(x$f_time_norm != 1)) {
+    cat("Time normalisation factors $f_time_norm:\n")
+    print(x$f_time_norm)
+  }
+  if (!is.null(x$meta)) {
+    cat("Meta information $meta:\n")
+    print(x$meta)
+  }
+  if (verbose) {
+    cat("\nDatasets $ds:")
+    for (ds in x$ds) {
+      cat("\n")
+      print(ds, data = data)
+    }
+  }
 }
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index f52692ba..388060d4 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean"))
+utils::globalVariables(c("name", "value_mean"))
 
 #' Calculate the minimum error to assume in order to pass the variance test
 #'
diff --git a/R/mkinerrplot.R b/R/mkinerrplot.R
index 36e22a43..02a12672 100644
--- a/R/mkinerrplot.R
+++ b/R/mkinerrplot.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual"))
+utils::globalVariables(c("variable", "residual"))
 
 #' Function to plot squared residuals and the error model for an mkin object
 #' 
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 7fa1c56e..a6efc858 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
+utils::globalVariables(c("name", "time", "value"))
 
 #' Fit a kinetic model to data with one or more state variables
 #'
diff --git a/R/mkinresplot.R b/R/mkinresplot.R
index 0bfdd02f..bad28ae8 100644
--- a/R/mkinresplot.R
+++ b/R/mkinresplot.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual"))
+utils::globalVariables(c("variable", "residual"))
 
 #' Function to plot residuals stored in an mkin object
 #'
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 903c3213..0d6b15c0 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables("ds")
+utils::globalVariables("ds")
 
 #' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
 #'
diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R
index 48df9483..6ac4e45d 100644
--- a/R/plot.mkinfit.R
+++ b/R/plot.mkinfit.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "observed"))
+utils::globalVariables(c("type", "variable", "observed"))
 
 #' Plot the observed data and the fitted model of an mkinfit object
 #'
diff --git a/R/saem.R b/R/saem.R
index f6d487ec..377da903 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -102,7 +102,7 @@ saem.mmkin <- function(object,
     grDevices::png(tmp)
   }
   fit_time <- system.time({
-    capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
+    utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
     f_pred <- try(saemix::saemix.predict(f_saemix), silent = TRUE)
     if (!inherits(f_pred, "try-error")) {
       f_saemix <- f_pred