aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/create_deg_func.R54
-rw-r--r--R/logistic.solution.R59
-rw-r--r--R/mkinfit.R2
-rw-r--r--R/mkinmod.R43
-rw-r--r--R/mkinpredict.R22
-rw-r--r--R/nlme.R6
-rw-r--r--R/parent_solutions.R4
7 files changed, 74 insertions, 116 deletions
diff --git a/R/create_deg_func.R b/R/create_deg_func.R
new file mode 100644
index 00000000..40efb3a3
--- /dev/null
+++ b/R/create_deg_func.R
@@ -0,0 +1,54 @@
+#' Create degradation functions for known analytical solutions
+#'
+#' @param spec List of model specifications as contained in mkinmod objects
+#' @param use_of_ff Minimum or maximum use of formation fractions
+#' @return Degradation function to be attached to mkinmod objects
+#' @examples
+#'
+#' SFO_SFO <- mkinmod(
+#' parent = mkinsub("SFO", "m1"),
+#' m1 = mkinsub("SFO"))
+#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
+
+create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
+
+ use_of_ff <- match.arg(use_of_ff)
+
+ min_ff <- use_of_ff == "min"
+
+ obs_vars <- names(spec)
+
+ n <- character(0)
+
+ parent <- obs_vars[1]
+
+ n[1] <- paste0(parent, " = ", spec[[1]]$type, ".solution(outtimes, odeini['", parent,
+ if (spec[[1]]$type == "SFORB") "_free", "'], ",
+ switch(spec[[1]]$type,
+ SFO = paste0("k_", parent, if (min_ff) "_sink" else "", ")"),
+ FOMC = "alpha, beta)",
+ IORE = paste0("k__iore_", parent, if (min_ff) "_sink" else "", ", N_", parent, ")"),
+ DFOP = "k1, k2, g)",
+ SFORB = paste0("k_", parent, "_free_bound, k_", parent, "_bound_free, k_", parent, "_free", if (min_ff) "_sink" else "", ")"),
+ HS = "k1, k2, tb)",
+ logistic = "kmax, k0, r)"
+ )
+ )
+
+ all_n <- paste(n, collapse = ",\n")
+
+ f_body <- paste0("{\n",
+ "out <- with(as.list(odeparms), {\n",
+ "data.frame(\n",
+ "time = outtimes,\n",
+ all_n, "\n",
+ ")})\n",
+ "return(out)\n}\n"
+ )
+
+ deg_func <- function(odeini, odeparms, outtimes) {}
+
+ body(deg_func) <- parse(text = f_body)
+
+ return(deg_func)
+}
diff --git a/R/logistic.solution.R b/R/logistic.solution.R
deleted file mode 100644
index d9db13d7..00000000
--- a/R/logistic.solution.R
+++ /dev/null
@@ -1,59 +0,0 @@
-#' Logistic kinetics
-#'
-#' Function describing exponential decline from a defined starting value, with
-#' an increasing rate constant, supposedly caused by microbial growth
-#'
-#' @param t Time.
-#' @param parent.0 Starting value for the response variable at time zero.
-#' @param kmax Maximum rate constant.
-#' @param k0 Minumum rate constant effective at time zero.
-#' @param r Growth rate of the increase in the rate constant.
-#' @return The value of the response variable at time \code{t}.
-#' @note The solution of the logistic model reduces to the
-#' \code{\link{SFO.solution}} if \code{k0} is equal to \code{kmax}.
-#' @references 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}
-#' @examples
-#'
-#' # Reproduce the plot on page 57 of FOCUS (2014)
-#' plot(function(x) logistic.solution(x, 100, 0.08, 0.0001, 0.2),
-#' from = 0, to = 100, ylim = c(0, 100),
-#' xlab = "Time", ylab = "Residue")
-#' plot(function(x) logistic.solution(x, 100, 0.08, 0.0001, 0.4),
-#' from = 0, to = 100, add = TRUE, lty = 2, col = 2)
-#' plot(function(x) logistic.solution(x, 100, 0.08, 0.0001, 0.8),
-#' from = 0, to = 100, add = TRUE, lty = 3, col = 3)
-#' plot(function(x) logistic.solution(x, 100, 0.08, 0.001, 0.2),
-#' from = 0, to = 100, add = TRUE, lty = 4, col = 4)
-#' plot(function(x) logistic.solution(x, 100, 0.08, 0.08, 0.2),
-#' from = 0, to = 100, add = TRUE, lty = 5, col = 5)
-#' legend("topright", inset = 0.05,
-#' legend = paste0("k0 = ", c(0.0001, 0.0001, 0.0001, 0.001, 0.08),
-#' ", r = ", c(0.2, 0.4, 0.8, 0.2, 0.2)),
-#' lty = 1:5, col = 1:5)
-#'
-#' # Fit with synthetic data
-#' logistic <- mkinmod(parent = mkinsub("logistic"))
-#'
-#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
-#' parms_logistic <- c(kmax = 0.08, k0 = 0.0001, r = 0.2)
-#' d_logistic <- mkinpredict(logistic,
-#' parms_logistic, c(parent = 100),
-#' sampling_times)
-#' d_2_1 <- add_err(d_logistic,
-#' sdfunc = function(x) sigma_twocomp(x, 0.5, 0.07),
-#' n = 1, reps = 2, digits = 5, LOD = 0.1, seed = 123456)[[1]]
-#'
-#' m <- mkinfit("logistic", d_2_1, quiet = TRUE)
-#' plot_sep(m)
-#' summary(m)$bpar
-#' endpoints(m)$distimes
-#'
-#' @export
-logistic.solution <- function(t, parent.0, kmax, k0, r)
-{
- parent = parent.0 * (kmax / (kmax - k0 + k0 * exp (r * t))) ^(kmax/r)
-}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 2b7e71cb..54dd75c2 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -269,7 +269,7 @@ mkinfit <- function(mkinmod, observed,
if (mkinmod[[1]] %in% parent_models_available) {
speclist <- list(list(type = mkinmod, sink = TRUE))
names(speclist) <- presumed_parent_name
- mkinmod <- mkinmod(speclist = speclist)
+ mkinmod <- mkinmod(speclist = speclist, use_of_ff = "min")
} else {
stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ",
paste(parent_models_available, collapse = ", "))
diff --git a/R/mkinmod.R b/R/mkinmod.R
index f52baa4f..21551861 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -101,7 +101,7 @@
#' }
#'
#' @export mkinmod
-mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verbose = FALSE)
+mkinmod <- function(..., use_of_ff = "max", speclist = NULL, quiet = FALSE, verbose = FALSE)
{
if (is.null(speclist)) spec <- list(...)
else spec <- speclist
@@ -421,45 +421,8 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb
}
# }}}
- # Attach a degradation function if an analytical solution is available {{{
- parent_type = spec[[1]]$type
- parent_name = names(spec)[[1]]
- if (length(spec) == 1) {
- odeparm_map <- switch(parent_type,
- SFO = c(
- k = if(use_of_ff == "min") paste("k", parent_name, "sink", sep = "_")
- else paste("k", parent_name, sep = "_")),
- FOMC = c(alpha = "alpha", beta = "beta"),
- IORE = c(
- k__iore = if(use_of_ff == "min") paste("k__iore", parent_name, "sink", sep = "_")
- else paste("k__iore", parent_name, sep = "_"),
- N = paste("N", parent_name, sep = "_")),
- DFOP = c(k1 = "k1", k2 = "k2", g = "g"),
- HS = c(k1 = "k1", k2 = "k2", tb = "tb"),
- SFORB = c(
- k_12 = paste("k", parent_name, "free_bound", sep = "_"),
- k_21 = paste("k", parent_name, "bound_free", sep = "_"),
- k_1output = paste("k", parent_name, "free_sink", sep = "_")),
- logistic = c(kmax = "kmax", k0 = "k0", r = "r")
- )
- odeparm_rev_map <- names(odeparm_map)
- names(odeparm_rev_map) <- odeparm_map
-
- model$deg_func <- function(odeini, odeparms, outtimes) {
- parent_func <- getFromNamespace(paste0(parent_type, ".solution"), "mkin")
- odeparm_list <- as.list(odeparms)
- names(odeparm_list) <- odeparm_rev_map[names(odeparm_list)]
-
- values <- do.call(parent_func,
- args = c(
- list(t = outtimes, parent_0 = odeini[1]),
- odeparm_list))
- out <- data.frame(outtimes, values)
- names(out) <- c("time", parent_name)
- return(out)
- }
- }
- # }}}
+ # Attach a degradation function if an analytical solution is available
+ model$deg_func <- create_deg_func(spec, use_of_ff)
class(model) <- "mkinmod"
return(model)
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 0f8e83bb..df51dbe3 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -43,36 +43,36 @@
#'
#' SFO <- mkinmod(degradinol = mkinsub("SFO"))
#' # Compare solution types
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' solution_type = "analytical")
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' solution_type = "deSolve")
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' solution_type = "deSolve", use_compiled = FALSE)
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' solution_type = "eigen")
#'
#' # Compare integration methods to analytical solution
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' solution_type = "analytical")[21,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' method = "lsoda")[21,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' method = "ode45")[21,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20,
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100), 0:20,
#' method = "rk4")[21,]
#' # rk4 is not as precise here
#'
#' # The number of output times used to make a lot of difference until the
#' # default for atol was adjusted
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100),
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100),
#' seq(0, 20, by = 0.1))[201,]
-#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100),
+#' mkinpredict(SFO, c(k_degradinol = 0.3), c(degradinol = 100),
#' seq(0, 20, by = 0.01))[2001,]
#'
#' # Check compiled model versions - they are faster than the eigenvalue based solutions!
#' SFO_SFO = mkinmod(parent = list(type = "SFO", to = "m1"),
-#' m1 = list(type = "SFO"))
+#' m1 = list(type = "SFO"), use_of_ff = "min")
#' if(require(rbenchmark)) {
#' benchmark(
#' eigen = mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01),
diff --git a/R/nlme.R b/R/nlme.R
index 60cd2ca1..fa69db3c 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -13,15 +13,15 @@
#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
#' m_SFO <- mkinmod(parent = mkinsub("SFO"))
#' d_SFO_1 <- mkinpredict(m_SFO,
-#' c(k_parent_sink = 0.1),
+#' c(k_parent = 0.1),
#' c(parent = 98), sampling_times)
#' d_SFO_1_long <- mkin_wide_to_long(d_SFO_1, time = "time")
#' d_SFO_2 <- mkinpredict(m_SFO,
-#' c(k_parent_sink = 0.05),
+#' c(k_parent = 0.05),
#' c(parent = 102), sampling_times)
#' d_SFO_2_long <- mkin_wide_to_long(d_SFO_2, time = "time")
#' d_SFO_3 <- mkinpredict(m_SFO,
-#' c(k_parent_sink = 0.02),
+#' c(k_parent = 0.02),
#' c(parent = 103), sampling_times)
#' d_SFO_3_long <- mkin_wide_to_long(d_SFO_3, time = "time")
#'
diff --git a/R/parent_solutions.R b/R/parent_solutions.R
index c33d6d13..e02bcda7 100644
--- a/R/parent_solutions.R
+++ b/R/parent_solutions.R
@@ -136,7 +136,7 @@ DFOP.solution <- function(t, parent_0, k1, k2, g)
#' between them.
#'
#' @family parent solutions
-#' @inherit HS.solution
+#' @inherit DFOP.solution
#' @param tb Break point. Before this time, exponential decline according to
#' \code{k1} is calculated, after this time, exponential decline proceeds
#' according to \code{k2}.
@@ -161,7 +161,7 @@ HS.solution <- function(t, parent_0, k1, k2, tb)
#' and no substance in the bound fraction.
#'
#' @family parent solutions
-#' @inherit HS.solution
+#' @inherit SFO.solution
#' @param k_12 Kinetic constant describing transfer from free to bound.
#' @param k_21 Kinetic constant describing transfer from bound to free.
#' @param k_1output Kinetic constant describing degradation of the free

Contact - Imprint