diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/create_deg_func.R | 54 | ||||
| -rw-r--r-- | R/logistic.solution.R | 59 | ||||
| -rw-r--r-- | R/mkinfit.R | 2 | ||||
| -rw-r--r-- | R/mkinmod.R | 43 | ||||
| -rw-r--r-- | R/mkinpredict.R | 22 | ||||
| -rw-r--r-- | R/nlme.R | 6 | ||||
| -rw-r--r-- | R/parent_solutions.R | 4 |
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), @@ -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 |