diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/AIC.mmkin.R | 50 | ||||
| -rw-r--r-- | R/CAKE_export.R | 44 | ||||
| -rw-r--r-- | R/DFOP.solution.R | 32 | ||||
| -rw-r--r-- | R/FOMC.solution.R | 40 | ||||
| -rw-r--r-- | R/HS.solution.R | 35 | ||||
| -rw-r--r-- | R/IORE.solution.R | 41 | ||||
| -rw-r--r-- | R/SFO.solution.R | 26 | ||||
| -rw-r--r-- | R/SFORB.solution.R | 44 | ||||
| -rw-r--r-- | R/add_err.R | 96 | ||||
| -rw-r--r-- | R/endpoints.R | 377 | ||||
| -rw-r--r-- | R/ilr.R | 46 | ||||
| -rw-r--r-- | R/logLik.mkinfit.R | 52 | ||||
| -rw-r--r-- | R/logistic.solution.R | 55 | ||||
| -rw-r--r-- | R/max_twa_parent.R | 70 | ||||
| -rw-r--r-- | R/mkin_long_to_wide.R | 57 | ||||
| -rw-r--r-- | R/mkin_wide_to_long.R | 68 | ||||
| -rw-r--r-- | R/mkinds.R | 55 | ||||
| -rw-r--r-- | R/mkinerrmin.R | 55 | ||||
| -rw-r--r-- | R/mkinerrplot.R | 56 | ||||
| -rw-r--r-- | R/mkinfit.R | 1801 | ||||
| -rw-r--r-- | R/mkinmod.R | 866 | ||||
| -rw-r--r-- | R/mkinparplot.R | 37 | ||||
| -rw-r--r-- | R/mkinplot.R | 4 | ||||
| -rw-r--r-- | R/mkinpredict.R | 119 | ||||
| -rw-r--r-- | R/mkinresplot.R | 149 | ||||
| -rw-r--r-- | R/mkinsub.R | 54 | ||||
| -rw-r--r-- | R/mmkin.R | 113 | ||||
| -rw-r--r-- | R/nafta.R | 79 | ||||
| -rw-r--r-- | R/plot.mkinfit.R | 126 | ||||
| -rw-r--r-- | R/plot.mmkin.R | 68 | ||||
| -rw-r--r-- | R/sigma_twocomp.R | 26 | ||||
| -rw-r--r-- | R/summary.mkinfit.R | 272 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 438 |
33 files changed, 3364 insertions, 2087 deletions
diff --git a/R/AIC.mmkin.R b/R/AIC.mmkin.R index 1d306ff9..ab17f0a0 100644 --- a/R/AIC.mmkin.R +++ b/R/AIC.mmkin.R @@ -1,21 +1,35 @@ -# Copyright (C) 2018 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/> -AIC.mmkin <- function(object, ..., k = 2) { +#' Calculated the AIC for a column of an mmkin object +#' +#' Provides a convenient way to compare different kinetic models fitted to the +#' same dataset. +#' +#' @param object An object of class \code{\link{mmkin}}, containing only one +#' column. +#' @param \dots For compatibility with the generic method +#' @param k As in the generic method +#' @return As in the generic method (a numeric value for single fits, or a +#' dataframe if there are several fits in the column). +#' @author Johannes Ranke +#' @examples +#' +#' \dontrun{ # skip, as it takes > 10 s on winbuilder +#' f <- mmkin(c("SFO", "FOMC", "DFOP"), +#' list("FOCUS A" = FOCUS_2006_A, +#' "FOCUS C" = FOCUS_2006_C), cores = 1, quiet = TRUE) +#' AIC(f[1, "FOCUS A"]) # We get a single number for a single fit +#' +#' # For FOCUS A, the models fit almost equally well, so the higher the number +#' # of parameters, the higher (worse) the AIC +#' AIC(f[, "FOCUS A"]) +#' AIC(f[, "FOCUS A"], k = 0) # If we do not penalize additional parameters, we get nearly the same +#' +#' # For FOCUS C, the more complex models fit better +#' AIC(f[, "FOCUS C"]) +#' } +#' +#' @export +AIC.mmkin <- function(object, ..., k = 2) +{ # We can only handle a single column if (ncol(object) != 1) stop("Please provide a single column object") n.fits <- length(object) diff --git a/R/CAKE_export.R b/R/CAKE_export.R index db9caa8d..70661b10 100644 --- a/R/CAKE_export.R +++ b/R/CAKE_export.R @@ -1,20 +1,30 @@ -# Copyright (C) 2019 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/> +#' Export a list of datasets format to a CAKE study file +#' +#' In addition to the datasets, the pathways in the degradation model can be +#' specified as well. +#' +#' @param ds A named list of datasets in long format as compatible with +#' \code{\link{mkinfit}}. +#' @param map A character vector with CAKE compartment names (Parent, A1, ...), +#' named with the names used in the list of datasets. +#' @param links An optional character vector of target compartments, named with +#' the names of the source compartments. In order to make this easier, the +#' names are used as in the datasets supplied. +#' @param filename Where to write the result. Should end in .csf in order to be +#' compatible with CAKE. +#' @param path An optional path to the output file. +#' @param overwrite If TRUE, existing files are overwritten. +#' @param study The name of the study. +#' @param description An optional description. +#' @param time_unit The time unit for the residue data. +#' @param res_unit The unit used for the residues. +#' @param comment An optional comment. +#' @param date The date of file creation. +#' @param optimiser Can be OLS or IRLS. +#' @importFrom utils write.table +#' @return The function is called for its side effect. +#' @author Johannes Ranke +#' @export CAKE_export <- function(ds, map = c(parent = "Parent"), links = NA, filename = "CAKE_export.csf", path = ".", overwrite = FALSE, diff --git a/R/DFOP.solution.R b/R/DFOP.solution.R index 8531cfed..608e7e18 100644 --- a/R/DFOP.solution.R +++ b/R/DFOP.solution.R @@ -1,5 +1,27 @@ -DFOP.solution <- function(t, parent.0, k1, k2, g)
-{
- parent = g * parent.0 * exp(-k1 * t) +
- (1 - g) * parent.0 * exp(-k2 * t)
-}
+#' Double First-Order in Parallel kinetics +#' +#' Function describing decline from a defined starting value using the sum of +#' two exponential decline functions. +#' +#' @param t Time. +#' @param parent.0 Starting value for the response variable at time zero. +#' @param k1 First kinetic constant. +#' @param k2 Second kinetic constant. +#' @param g Fraction of the starting value declining according to the first +#' kinetic constant. +#' @return The value of the response variable at time \code{t}. +#' @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} +#' @examples +#' +#' plot(function(x) DFOP.solution(x, 100, 5, 0.5, 0.3), 0, 4, ylim = c(0,100)) +#' +#' @export +DFOP.solution <- function(t, parent.0, k1, k2, g) +{ + parent = g * parent.0 * exp(-k1 * t) + + (1 - g) * parent.0 * exp(-k2 * t) +} diff --git a/R/FOMC.solution.R b/R/FOMC.solution.R index 8bd13d6b..f5e6a7ea 100644 --- a/R/FOMC.solution.R +++ b/R/FOMC.solution.R @@ -1,4 +1,36 @@ -FOMC.solution <- function(t, parent.0, alpha, beta)
-{
- parent = parent.0 / (t/beta + 1)^alpha
-}
+#' First-Order Multi-Compartment kinetics +#' +#' Function describing exponential decline from a defined starting value, with +#' a decreasing rate constant. +#' +#' The form given here differs slightly from the original reference by +#' Gustafson and Holden (1990). The parameter \code{beta} corresponds to 1/beta +#' in the original equation. +#' +#' @param t Time. +#' @param parent.0 Starting value for the response variable at time zero. +#' @param alpha Shape parameter determined by coefficient of variation of rate +#' constant values. +#' @param beta Location parameter. +#' @return The value of the response variable at time \code{t}. +#' @note The solution of the FOMC kinetic model reduces to the +#' \code{\link{SFO.solution}} for large values of \code{alpha} and +#' \code{beta} with \eqn{k = \frac{\beta}{\alpha}}{k = beta/alpha}. +#' @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} +#' +#' Gustafson DI and Holden LR (1990) Nonlinear pesticide dissipation in soil: +#' A new model based on spatial variability. \emph{Environmental Science and +#' Technology} \bold{24}, 1032-1038 +#' @examples +#' +#' plot(function(x) FOMC.solution(x, 100, 10, 2), 0, 2, ylim = c(0, 100)) +#' +#' @export +FOMC.solution <- function(t, parent.0, alpha, beta) +{ + parent = parent.0 / (t/beta + 1)^alpha +} diff --git a/R/HS.solution.R b/R/HS.solution.R index 4651a6a8..890ad8ff 100644 --- a/R/HS.solution.R +++ b/R/HS.solution.R @@ -1,6 +1,29 @@ -HS.solution <- function(t, parent.0, k1, k2, tb)
-{
- parent = ifelse(t <= tb,
- parent.0 * exp(-k1 * t),
- parent.0 * exp(-k1 * tb) * exp(-k2 * (t - tb)))
-}
+#' Hockey-Stick kinetics +#' +#' Function describing two exponential decline functions with a break point +#' between them. +#' +#' @param t Time. +#' @param parent.0 Starting value for the response variable at time zero. +#' @param k1 First kinetic constant. +#' @param k2 Second kinetic constant. +#' @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}. +#' @return The value of the response variable at time \code{t}. +#' @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} +#' @examples +#' +#' plot(function(x) HS.solution(x, 100, 2, 0.3, 0.5), 0, 2, ylim=c(0,100)) +#' +#' @export +HS.solution <- function(t, parent.0, k1, k2, tb) +{ + parent = ifelse(t <= tb, + parent.0 * exp(-k1 * t), + parent.0 * exp(-k1 * tb) * exp(-k2 * (t - tb))) +} diff --git a/R/IORE.solution.R b/R/IORE.solution.R index 5405be96..58807108 100644 --- a/R/IORE.solution.R +++ b/R/IORE.solution.R @@ -1,4 +1,37 @@ -IORE.solution <- function(t, parent.0, k__iore, N)
-{
- parent = (parent.0^(1 - N) - (1 - N) * k__iore * t)^(1/(1 - N))
-}
+#' Indeterminate order rate equation kinetics +#' +#' Function describing exponential decline from a defined starting value, with +#' a concentration dependent rate constant. +#' +#' @param t Time. +#' @param parent.0 Starting value for the response variable at time zero. +#' @param k__iore Rate constant. Note that this depends on the concentration +#' units used. +#' @param N Exponent describing the nonlinearity of the rate equation +#' @return The value of the response variable at time \code{t}. +#' @note The solution of the IORE kinetic model reduces to the +#' \code{\link{SFO.solution}} if N = 1. The parameters of the IORE model can +#' be transformed to equivalent parameters of the FOMC mode - see the NAFTA +#' guidance for details. +#' @references NAFTA Technical Working Group on Pesticides (not dated) Guidance +#' for Evaluating and Calculating Degradation Kinetics in Environmental Media +#' @keywords manip +#' @examples +#' +#' plot(function(x) IORE.solution(x, 100, 0.2, 1.3), 0, 2, ylim = c(0, 100)) +#' \dontrun{ +#' fit.fomc <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) +#' fit.iore <- mkinfit("IORE", FOCUS_2006_C, quiet = TRUE) +#' fit.iore.deS <- mkinfit("IORE", FOCUS_2006_C, solution_type = "deSolve", quiet = TRUE) +#' +#' print(data.frame(fit.fomc$par, fit.iore$par, fit.iore.deS$par, +#' row.names = paste("model par", 1:4))) +#' print(rbind(fomc = endpoints(fit.fomc)$distimes, iore = endpoints(fit.iore)$distimes, +#' iore.deS = endpoints(fit.iore)$distimes)) +#' } +#' +#' @export +IORE.solution <- function(t, parent.0, k__iore, N) +{ + parent = (parent.0^(1 - N) - (1 - N) * k__iore * t)^(1/(1 - N)) +} diff --git a/R/SFO.solution.R b/R/SFO.solution.R index 3a376e48..17c16a4d 100644 --- a/R/SFO.solution.R +++ b/R/SFO.solution.R @@ -1,4 +1,22 @@ -SFO.solution <- function(t, parent.0, k)
-{
- parent = parent.0 * exp(-k * t)
-}
+#' Single First-Order kinetics +#' +#' Function describing exponential decline from a defined starting value. +#' +#' @param t Time. +#' @param parent.0 Starting value for the response variable at time zero. +#' @param k Kinetic constant. +#' @return The value of the response variable at time \code{t}. +#' @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} +#' @examples +#' +#' \dontrun{plot(function(x) SFO.solution(x, 100, 3), 0, 2)} +#' +#' @export +SFO.solution <- function(t, parent.0, k) +{ + parent = parent.0 * exp(-k * t) +} diff --git a/R/SFORB.solution.R b/R/SFORB.solution.R index 4cb94def..2abe4577 100644 --- a/R/SFORB.solution.R +++ b/R/SFORB.solution.R @@ -1,9 +1,35 @@ -SFORB.solution = function(t, parent.0, k_12, k_21, k_1output) {
- sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
- b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
- b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
-
- parent = parent.0 *
- (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
- ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
-}
+#' Single First-Order Reversible Binding kinetics +#' +#' Function describing the solution of the differential equations describing +#' the kinetic model with first-order terms for a two-way transfer from a free +#' to a bound fraction, and a first-order degradation term for the free +#' fraction. The initial condition is a defined amount in the free fraction +#' and no substance in the bound fraction. +#' +#' @param t Time. +#' @param parent.0 Starting value for the response variable at time zero. +#' @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 +#' fraction. +#' @return The value of the response variable, which is the sum of free and +#' bound fractions at time \code{t}. +#' @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} +#' @examples +#' +#' \dontrun{plot(function(x) SFORB.solution(x, 100, 0.5, 2, 3), 0, 2)} +#' +#' @export +SFORB.solution = function(t, parent.0, k_12, k_21, k_1output) { + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + parent = parent.0 * + (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) +} diff --git a/R/add_err.R b/R/add_err.R index b2a1808e..a523e9c2 100644 --- a/R/add_err.R +++ b/R/add_err.R @@ -1,24 +1,78 @@ -# Copyright (C) 2015-2018 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/> - -add_err = function(prediction, sdfunc, secondary = c("M1", "M2"), - n = 1000, LOD = 0.1, reps = 2, - digits = 1, seed = NA) +#' Add normally distributed errors to simulated kinetic degradation data +#' +#' Normally distributed errors are added to data predicted for a specific +#' degradation model using \code{\link{mkinpredict}}. The variance of the error +#' may depend on the predicted value and is specified as a standard deviation. +#' +#' @param prediction A prediction from a kinetic model as produced by +#' \code{\link{mkinpredict}}. +#' @param sdfunc A function taking the predicted value as its only argument and +#' returning a standard deviation that should be used for generating the +#' random error terms for this value. +#' @param secondary The names of state variables that should have an initial +#' value of zero +#' @param n The number of datasets to be generated. +#' @param LOD The limit of detection (LOD). Values that are below the LOD after +#' adding the random error will be set to NA. +#' @param reps The number of replicates to be generated within the datasets. +#' @param digits The number of digits to which the values will be rounded. +#' @param seed The seed used for the generation of random numbers. If NA, the +#' seed is not set. +#' @importFrom stats rnorm +#' @return A list of datasets compatible with \code{\link{mmkin}}, i.e. the +#' components of the list are datasets compatible with \code{\link{mkinfit}}. +#' @author Johannes Ranke +#' @references Ranke J and Lehmann R (2015) To t-test or not to t-test, that is +#' the question. XV Symposium on Pesticide Chemistry 2-4 September 2015, +#' Piacenza, Italy +#' http://chem.uft.uni-bremen.de/ranke/posters/piacenza_2015.pdf +#' @examples +#' +#' # The kinetic model +#' m_SFO_SFO <- mkinmod(parent = mkinsub("SFO", "M1"), +#' M1 = mkinsub("SFO"), use_of_ff = "max") +#' +#' # Generate a prediction for a specific set of parameters +#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +#' +#' # This is the prediction used for the "Type 2 datasets" on the Piacenza poster +#' # from 2015 +#' d_SFO_SFO <- mkinpredict(m_SFO_SFO, +#' c(k_parent = 0.1, f_parent_to_M1 = 0.5, +#' k_M1 = log(2)/1000), +#' c(parent = 100, M1 = 0), +#' sampling_times) +#' +#' # Add an error term with a constant (independent of the value) standard deviation +#' # of 10, and generate three datasets +#' d_SFO_SFO_err <- add_err(d_SFO_SFO, function(x) 10, n = 3, seed = 123456789 ) +#' +#' # Name the datasets for nicer plotting +#' names(d_SFO_SFO_err) <- paste("Dataset", 1:3) +#' +#' # Name the model in the list of models (with only one member in this case) for +#' # nicer plotting later on. Be quiet and use only one core not to offend CRAN +#' # checks +#' \dontrun{ +#' f_SFO_SFO <- mmkin(list("SFO-SFO" = m_SFO_SFO), +#' d_SFO_SFO_err, cores = 1, +#' quiet = TRUE) +#' +#' plot(f_SFO_SFO) +#' +#' # We would like to inspect the fit for dataset 3 more closely +#' # Using double brackets makes the returned object an mkinfit object +#' # instead of a list of mkinfit objects, so plot.mkinfit is used +#' plot(f_SFO_SFO[[3]], show_residuals = TRUE) +#' +#' # If we use single brackets, we should give two indices (model and dataset), +#' # and plot.mmkin is used +#' plot(f_SFO_SFO[1, 3]) +#' } +#' +#' @export +add_err <- function(prediction, sdfunc, secondary = c("M1", "M2"), + n = 1000, LOD = 0.1, reps = 2, digits = 1, seed = NA) { if (!is.na(seed)) set.seed(seed) diff --git a/R/endpoints.R b/R/endpoints.R index f8a44c4d..14beadea 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -1,178 +1,199 @@ -endpoints <- function(fit) {
- # Calculate dissipation times DT50 and DT90 and formation
- # fractions as well as SFORB eigenvalues from optimised parameters
- # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
- # HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
- ep <- list()
- obs_vars <- fit$obs_vars
- parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
- ep$ff <- vector()
- ep$SFORB <- vector()
- ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
- DT90 = rep(NA, length(obs_vars)),
- row.names = obs_vars)
- for (obs_var in obs_vars) {
- type = names(fit$mkinmod$map[[obs_var]])[1]
-
- # Get formation fractions if directly fitted, and calculate remaining fraction to sink
- f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE)
- if (length(f_names) > 0) {
- f_values = parms.all[f_names]
- f_to_sink = 1 - sum(f_values)
- names(f_to_sink) = ifelse(type == "SFORB",
- paste(obs_var, "free", "sink", sep = "_"),
- paste(obs_var, "sink", sep = "_"))
- for (f_name in f_names) {
- ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
- }
- ep$ff = append(ep$ff, f_to_sink)
- }
-
- # Get the rest
- if (type == "SFO") {
- k_names = grep(paste("^k", obs_var, sep="_"), names(parms.all), value=TRUE)
- k_tot = sum(parms.all[k_names])
- DT50 = log(2)/k_tot
- DT90 = log(10)/k_tot
- if (fit$mkinmod$use_of_ff == "min") {
- for (k_name in k_names)
- {
- ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
- }
- }
- }
- if (type == "FOMC") {
- alpha = parms.all["alpha"]
- beta = parms.all["beta"]
- DT50 = beta * (2^(1/alpha) - 1)
- DT90 = beta * (10^(1/alpha) - 1)
- DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
- ep$distimes[obs_var, c("DT50back")] = DT50_back
- }
- if (type == "IORE") {
- k_names = grep(paste("^k__iore", obs_var, sep="_"), names(parms.all), value=TRUE)
- k_tot = sum(parms.all[k_names])
- # From the NAFTA kinetics guidance, p. 5
- n = parms.all[paste("N", obs_var, sep = "_")]
- k = k_tot
- # Use the initial concentration of the parent compound
- source_name = fit$mkinmod$map[[1]][[1]]
- c0 = parms.all[paste(source_name, "0", sep = "_")]
- alpha = 1 / (n - 1)
- beta = (c0^(1 - n))/(k * (n - 1))
- DT50 = beta * (2^(1/alpha) - 1)
- DT90 = beta * (10^(1/alpha) - 1)
- DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
- ep$distimes[obs_var, c("DT50back")] = DT50_back
- if (fit$mkinmod$use_of_ff == "min") {
- for (k_name in k_names)
- {
- ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
- }
- }
- }
- if (type == "DFOP") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- g = parms.all["g"]
- f <- function(log_t, x) {
- t <- exp(log_t)
- fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
- (fraction - (1 - x/100))^2
- }
- DT50_k1 = log(2)/k1
- DT50_k2 = log(2)/k2
- DT90_k1 = log(10)/k1
- DT90_k2 = log(10)/k2
-
- DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum),
- silent = TRUE)
- DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum),
- silent = TRUE)
- if (inherits(DT50, "try-error")) DT50 = NA
- if (inherits(DT90, "try-error")) DT90 = NA
-
- ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
- ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
- }
- if (type == "HS") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- tb = parms.all["tb"]
- DTx <- function(x) {
- DTx.a <- (log(100/(100 - x)))/k1
- DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
- if (DTx.a < tb) DTx <- DTx.a
- else DTx <- DTx.b
- return(DTx)
- }
- DT50 <- DTx(50)
- DT90 <- DTx(90)
- DT50_k1 = log(2)/k1
- DT50_k2 = log(2)/k2
- ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
- ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
- }
- if (type == "SFORB") {
- # FOCUS kinetics (2006), p. 60 f
- k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
- k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
- k_1output = sum(parms.all[k_out_names])
- k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]
- k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]
-
- sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
- b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
- b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
-
- DT50_b1 = log(2)/b1
- DT50_b2 = log(2)/b2
- DT90_b1 = log(10)/b1
- DT90_b2 = log(10)/b2
-
- SFORB_fraction = function(t) {
- ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
- ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
- }
-
- f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2
- log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum,
- silent = TRUE)
- f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2
- log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum,
- silent = TRUE)
-
- DT50 = if (inherits(log_DT50, "try-error")) NA
- else exp(log_DT50)
- DT90 = if (inherits(log_DT90, "try-error")) NA
- else exp(log_DT90)
-
- for (k_out_name in k_out_names)
- {
- ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output
- }
-
- # Return the eigenvalues for comparison with DFOP rate constants
- ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
- ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
-
- ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1
- ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2
- }
- if (type == "logistic") {
- # FOCUS kinetics (2014) p. 67
- kmax = parms.all["kmax"]
- k0 = parms.all["k0"]
- r = parms.all["r"]
- DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))
- DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))
-
- DT50_k0 = log(2)/k0
- DT50_kmax = log(2)/kmax
- ep$distimes[obs_var, c("DT50_k0")] = DT50_k0
- ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax
- }
- ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)
- }
- return(ep)
-}
+#' Function to calculate endpoints for further use from kinetic models fitted +#' with mkinfit +#' +#' This function calculates DT50 and DT90 values as well as formation fractions +#' from kinetic models fitted with mkinfit. If the SFORB model was specified +#' for one of the parents or metabolites, the Eigenvalues are returned. These +#' are equivalent to the rate constantes of the DFOP model, but with the +#' advantage that the SFORB model can also be used for metabolites. +#' +#' @param fit An object of class \code{\link{mkinfit}}. +#' @importFrom stats optimize +#' @return A list with the components mentioned above. +#' @note The function is used internally by \code{\link{summary.mkinfit}}. +#' @author Johannes Ranke +#' @keywords manip +#' @examples +#' +#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) +#' endpoints(fit) +#' +#' @export +endpoints <- function(fit) { + # Calculate dissipation times DT50 and DT90 and formation + # fractions as well as SFORB eigenvalues from optimised parameters + # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from + # HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models + ep <- list() + obs_vars <- fit$obs_vars + parms.all <- c(fit$bparms.optim, fit$bparms.fixed) + ep$ff <- vector() + ep$SFORB <- vector() + ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), + DT90 = rep(NA, length(obs_vars)), + row.names = obs_vars) + for (obs_var in obs_vars) { + type = names(fit$mkinmod$map[[obs_var]])[1] + + # Get formation fractions if directly fitted, and calculate remaining fraction to sink + f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE) + if (length(f_names) > 0) { + f_values = parms.all[f_names] + f_to_sink = 1 - sum(f_values) + names(f_to_sink) = ifelse(type == "SFORB", + paste(obs_var, "free", "sink", sep = "_"), + paste(obs_var, "sink", sep = "_")) + for (f_name in f_names) { + ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]] + } + ep$ff = append(ep$ff, f_to_sink) + } + + # Get the rest + if (type == "SFO") { + k_names = grep(paste("^k", obs_var, sep="_"), names(parms.all), value=TRUE) + k_tot = sum(parms.all[k_names]) + DT50 = log(2)/k_tot + DT90 = log(10)/k_tot + if (fit$mkinmod$use_of_ff == "min") { + for (k_name in k_names) + { + ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot + } + } + } + if (type == "FOMC") { + alpha = parms.all["alpha"] + beta = parms.all["beta"] + DT50 = beta * (2^(1/alpha) - 1) + DT90 = beta * (10^(1/alpha) - 1) + DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011 + ep$distimes[obs_var, c("DT50back")] = DT50_back + } + if (type == "IORE") { + k_names = grep(paste("^k__iore", obs_var, sep="_"), names(parms.all), value=TRUE) + k_tot = sum(parms.all[k_names]) + # From the NAFTA kinetics guidance, p. 5 + n = parms.all[paste("N", obs_var, sep = "_")] + k = k_tot + # Use the initial concentration of the parent compound + source_name = fit$mkinmod$map[[1]][[1]] + c0 = parms.all[paste(source_name, "0", sep = "_")] + alpha = 1 / (n - 1) + beta = (c0^(1 - n))/(k * (n - 1)) + DT50 = beta * (2^(1/alpha) - 1) + DT90 = beta * (10^(1/alpha) - 1) + DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011 + ep$distimes[obs_var, c("DT50back")] = DT50_back + if (fit$mkinmod$use_of_ff == "min") { + for (k_name in k_names) + { + ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot + } + } + } + if (type == "DFOP") { + k1 = parms.all["k1"] + k2 = parms.all["k2"] + g = parms.all["g"] + f <- function(log_t, x) { + t <- exp(log_t) + fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t) + (fraction - (1 - x/100))^2 + } + DT50_k1 = log(2)/k1 + DT50_k2 = log(2)/k2 + DT90_k1 = log(10)/k1 + DT90_k2 = log(10)/k2 + + DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum), + silent = TRUE) + DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum), + silent = TRUE) + if (inherits(DT50, "try-error")) DT50 = NA + if (inherits(DT90, "try-error")) DT90 = NA + + ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 + ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 + } + if (type == "HS") { + k1 = parms.all["k1"] + k2 = parms.all["k2"] + tb = parms.all["tb"] + DTx <- function(x) { + DTx.a <- (log(100/(100 - x)))/k1 + DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2 + if (DTx.a < tb) DTx <- DTx.a + else DTx <- DTx.b + return(DTx) + } + DT50 <- DTx(50) + DT90 <- DTx(90) + DT50_k1 = log(2)/k1 + DT50_k2 = log(2)/k2 + ep$distimes[obs_var, c("DT50_k1")] = DT50_k1 + ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 + } + if (type == "SFORB") { + # FOCUS kinetics (2006), p. 60 f + k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(parms.all), value=TRUE) + k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_")) + k_1output = sum(parms.all[k_out_names]) + k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")] + k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")] + + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + DT50_b1 = log(2)/b1 + DT50_b2 = log(2)/b2 + DT90_b1 = log(10)/b1 + DT90_b2 = log(10)/b2 + + SFORB_fraction = function(t) { + ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t) + } + + f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2 + log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum, + silent = TRUE) + f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2 + log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum, + silent = TRUE) + + DT50 = if (inherits(log_DT50, "try-error")) NA + else exp(log_DT50) + DT90 = if (inherits(log_DT90, "try-error")) NA + else exp(log_DT90) + + for (k_out_name in k_out_names) + { + ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output + } + + # Return the eigenvalues for comparison with DFOP rate constants + ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1 + ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2 + + ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1 + ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2 + } + if (type == "logistic") { + # FOCUS kinetics (2014) p. 67 + kmax = parms.all["kmax"] + k0 = parms.all["k0"] + r = parms.all["r"] + DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax)))) + DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax)))) + + DT50_k0 = log(2)/k0 + DT50_kmax = log(2)/kmax + ep$distimes[obs_var, c("DT50_k0")] = DT50_k0 + ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax + } + ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90) + } + return(ep) +} @@ -1,6 +1,3 @@ -# Copyright (C) 2012 René Lehmann, 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 @@ -16,6 +13,47 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see <http://www.gnu.org/licenses/> +#' Function to perform isometric log-ratio transformation +#' +#' This implementation is a special case of the class of isometric log-ratio +#' transformations. +#' +#' @aliases ilr invilr +#' @param x A numeric vector. Naturally, the forward transformation is only +#' sensible for vectors with all elements being greater than zero. +#' @return The result of the forward or backward transformation. The returned +#' components always sum to 1 for the case of the inverse log-ratio +#' transformation. +#' @author René Lehmann and Johannes Ranke +#' @seealso Another implementation can be found in R package +#' \code{robCompositions}. +#' @references Peter Filzmoser, Karel Hron (2008) Outlier Detection for +#' Compositional Data Using Robust Methods. Math Geosci 40 233-248 +#' @keywords manip +#' @examples +#' +#' # Order matters +#' ilr(c(0.1, 1, 10)) +#' ilr(c(10, 1, 0.1)) +#' # Equal entries give ilr transformations with zeros as elements +#' ilr(c(3, 3, 3)) +#' # Almost equal entries give small numbers +#' ilr(c(0.3, 0.4, 0.3)) +#' # Only the ratio between the numbers counts, not their sum +#' invilr(ilr(c(0.7, 0.29, 0.01))) +#' invilr(ilr(2.1 * c(0.7, 0.29, 0.01))) +#' # Inverse transformation of larger numbers gives unequal elements +#' invilr(-10) +#' invilr(c(-10, 0)) +#' # The sum of the elements of the inverse ilr is 1 +#' sum(invilr(c(-10, 0))) +#' # This is why we do not need all elements of the inverse transformation to go back: +#' a <- c(0.1, 0.3, 0.5) +#' b <- invilr(a) +#' length(b) # Four elements +#' ilr(c(b[1:3], 1 - sum(b[1:3]))) # Gives c(0.1, 0.3, 0.5) +#' +#' @export ilr <- function(x) { z <- vector() for (i in 1:(length(x) - 1)) { @@ -24,6 +62,8 @@ ilr <- function(x) { return(z) } +#' @rdname ilr +#' @export invilr<-function(x) { D <- length(x) + 1 z <- c(x, 0) diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R index d812f177..4ec3d7d4 100644 --- a/R/logLik.mkinfit.R +++ b/R/logLik.mkinfit.R @@ -1,24 +1,40 @@ -# Copyright (C) 2018,2019 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/> +#' Calculated the log-likelihood of a fitted mkinfit object +#' +#' This function simply calculates the product of the likelihood densities +#' calculated using \code{\link{dnorm}}, i.e. assuming normal distribution, +#' with of the mean predicted by the degradation model, and the standard +#' deviation predicted by the error model. +#' +#' The total number of estimated parameters returned with the value of the +#' likelihood is calculated as the sum of fitted degradation model parameters +#' and the fitted error model parameters. +#' +#' @param object An object of class \code{\link{mkinfit}}. +#' @param \dots For compatibility with the generic method +#' @return An object of class \code{\link{logLik}} with the number of estimated +#' parameters (degradation model parameters plus variance model parameters) +#' as attribute. +#' @author Johannes Ranke +#' @seealso Compare the AIC of columns of \code{\link{mmkin}} objects using +#' \code{\link{AIC.mmkin}}. +#' @examples +#' +#' \dontrun{ +#' sfo_sfo <- mkinmod( +#' parent = mkinsub("SFO", to = "m1"), +#' m1 = mkinsub("SFO") +#' ) +#' d_t <- FOCUS_2006_D +#' f_nw <- mkinfit(sfo_sfo, d_t, quiet = TRUE) # no weighting (weights are unity) +#' f_obs <- mkinfit(sfo_sfo, d_t, error_model = "obs", quiet = TRUE) +#' f_tc <- mkinfit(sfo_sfo, d_t, error_model = "tc", quiet = TRUE) +#' AIC(f_nw, f_obs, f_tc) +#' } +#' +#' @export logLik.mkinfit <- function(object, ...) { val <- object$logLik # Number of estimated parameters attr(val, "df") <- length(object$bparms.optim) + length(object$errparms) return(val) } -# vim: set ts=2 sw=2 expandtab: diff --git a/R/logistic.solution.R b/R/logistic.solution.R index a3bddab3..d9db13d7 100644 --- a/R/logistic.solution.R +++ b/R/logistic.solution.R @@ -1,3 +1,58 @@ +#' 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/max_twa_parent.R b/R/max_twa_parent.R index 5129e369..ef3c0ada 100644 --- a/R/max_twa_parent.R +++ b/R/max_twa_parent.R @@ -1,21 +1,43 @@ -# Copyright (C) 2016-2019 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/> - +#' Function to calculate maximum time weighted average concentrations from +#' kinetic models fitted with mkinfit +#' +#' 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, DFOP and HS models, using the analytical formulas given in the PEC +#' soil section of the FOCUS guidance. +#' +#' @aliases max_twa_parent max_twa_sfo max_twa_fomc max_twa_dfop max_twa_hs +#' @param fit An object of class \code{\link{mkinfit}}. +#' @param windows The width of the time windows for which the TWAs should be +#' calculated. +#' @param M0 The initial concentration for which the maximum time weighted +#' average over the decline curve should be calculated. The default is to use +#' a value of 1, which means that a relative maximum time weighted average +#' factor (f_twa) is calculated. +#' @param k The rate constant in the case of SFO kinetics. +#' @param t The width of the time window. +#' @param alpha Parameter of the FOMC model. +#' @param beta Parameter of the FOMC model. +#' @param k1 The first rate constant of the DFOP or the HS kinetics. +#' @param k2 The second rate constant of the DFOP or the HS kinetics. +#' @param g Parameter of the DFOP model. +#' @param tb Parameter of the HS model. +#' @return For \code{max_twa_parent}, a numeric vector, named using the +#' \code{windows} argument. For the other functions, a numeric vector of +#' length one (also known as 'a number'). +#' @author Johannes Ranke +#' @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} +#' @examples +#' +#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) +#' max_twa_parent(fit, c(7, 21)) +#' +#' @export max_twa_parent <- function(fit, windows) { parms.all <- c(fit$bparms.optim, fit$bparms.fixed) obs_vars <- fit$obs_vars @@ -74,15 +96,27 @@ max_twa_parent <- function(fit, windows) { names(res) <- windows return(res) } + +#' @rdname max_twa_parent +#' @export max_twa_sfo <- function(M0 = 1, k, t) { M0 * (1 - exp(- k * t)) / (k * t) } + +#' @rdname max_twa_parent +#' @export max_twa_fomc <- function(M0 = 1, alpha, beta, t) { M0 * (beta)/(t * (1 - alpha)) * ((t/beta + 1)^(1 - alpha) - 1) } + +#' @rdname max_twa_parent +#' @export max_twa_dfop <- function(M0 = 1, k1, k2, g, t) { M0/t * ((g/k1) * (1 - exp(- k1 * t)) + ((1 - g)/k2) * (1 - exp(- k2 * t))) } + +#' @rdname max_twa_parent +#' @export max_twa_hs <- function(M0 = 1, k1, k2, tb, t) { (M0 / t) * ( (1/k1) * (1 - exp(- k1 * tb)) + diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R index 081262f8..0125f2da 100644 --- a/R/mkin_long_to_wide.R +++ b/R/mkin_long_to_wide.R @@ -1,28 +1,29 @@ -# Copyright (C) 2010-2011 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/>
-
-mkin_long_to_wide <- function(long_data, time = "time", outtime = "time")
-{
- colnames <- unique(long_data$name)
- wide_data <- data.frame(time = subset(long_data, name == colnames[1], time))
- names(wide_data) <- outtime
- for (var in colnames) {
- wide_data[var] <- subset(long_data, name == var, value)
- }
- return(wide_data)
-}
+#' Convert a dataframe from long to wide format +#' +#' This function takes a dataframe in the long form, i.e. with a row for each +#' observed value, and converts it into a dataframe with one independent +#' variable and several dependent variables as columns. +#' +#' @param long_data The dataframe must contain one variable called "time" with +#' the time values specified by the \code{time} argument, one column called +#' "name" with the grouping of the observed values, and finally one column of +#' observed values called "value". +#' @param time The name of the time variable in the long input data. +#' @param outtime The name of the time variable in the wide output data. +#' @return Dataframe in wide format. +#' @author Johannes Ranke +#' @examples +#' +#' mkin_long_to_wide(FOCUS_2006_D) +#' +#' @export mkin_long_to_wide +mkin_long_to_wide <- function(long_data, time = "time", outtime = "time") +{ + colnames <- unique(long_data$name) + wide_data <- data.frame(time = subset(long_data, name == colnames[1], time)) + names(wide_data) <- outtime + for (var in colnames) { + wide_data[var] <- subset(long_data, name == var, value) + } + return(wide_data) +} diff --git a/R/mkin_wide_to_long.R b/R/mkin_wide_to_long.R index 21a7cbe9..bef0e408 100644 --- a/R/mkin_wide_to_long.R +++ b/R/mkin_wide_to_long.R @@ -1,34 +1,34 @@ -# $Id$
-
-# Copyright (C) 2010-2013 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
-
-mkin_wide_to_long <- function(wide_data, time = "t")
-{
- colnames <- names(wide_data)
- if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".")
- vars <- subset(colnames, colnames != time)
- n <- length(colnames) - 1
- long_data <- data.frame(
- name = rep(vars, each = length(wide_data[[time]])),
- time = as.numeric(rep(wide_data[[time]], n)),
- value = as.numeric(unlist(wide_data[vars])),
- row.names = NULL)
- return(long_data)
-}
+if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) + +#' Convert a dataframe with observations over time into long format +#' +#' This function simply takes a dataframe with one independent variable and +#' several dependent variable and converts it into the long form as required by +#' \code{\link{mkinfit}}. +#' +#' @param wide_data The dataframe must contain one variable with the time +#' values specified by the \code{time} argument and usually more than one +#' column of observed values. +#' @param time The name of the time variable. +#' @return Dataframe in long format as needed for \code{\link{mkinfit}}. +#' @author Johannes Ranke +#' @keywords manip +#' @examples +#' +#' wide <- data.frame(t = c(1,2,3), x = c(1,4,7), y = c(3,4,5)) +#' mkin_wide_to_long(wide) +#' +#' @export +mkin_wide_to_long <- function(wide_data, time = "t") +{ + colnames <- names(wide_data) + if (!(time %in% colnames)) stop("The data in wide format have to contain a variable named ", time, ".") + vars <- subset(colnames, colnames != time) + n <- length(colnames) - 1 + long_data <- data.frame( + name = rep(vars, each = length(wide_data[[time]])), + time = as.numeric(rep(wide_data[[time]], n)), + value = as.numeric(unlist(wide_data[vars])), + row.names = NULL) + return(long_data) +} @@ -1,21 +1,33 @@ -# Copyright (C) 2015,2018,2019 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/> - +#' A dataset class for mkin +#' +#' A dataset class for mkin +#' +#' @name mkinds +#' @docType class +#' @format An \code{\link{R6Class}} generator object. +#' @section Fields: +#' +#' \describe{ \item{list("title")}{A full title for the dataset} +#' +#' \item{list("sampling")}{times The sampling times} +#' +#' \item{list("time_unit")}{The time unit} +#' +#' \item{list("observed")}{Names of the observed compounds} +#' +#' \item{list("unit")}{The unit of the observations} +#' +#' \item{list("replicates")}{The number of replicates} +#' +#' \item{list("data")}{A dataframe with at least the columns name, time and +#' value in order to be compatible with mkinfit} } +#' @importFrom R6 R6Class +#' @keywords datasets +#' @examples +#' +#' mds <- mkinds$new("FOCUS A", FOCUS_2006_A) +#' +#' @export mkinds <- R6Class("mkinds", public = list( title = NULL, @@ -42,6 +54,13 @@ mkinds <- R6Class("mkinds", ) ) +#' Print mkinds objects +#' +#' Print mkinds objects. +#' +#' @param x An \code{\link{mkinds}} object. +#' @param \dots Not used. +#' @export print.mkinds <- function(x, ...) { cat("<mkinds> with $title: ", x$title, "\n") cat("Observed compounds $observed: ", paste(x$observed, collapse = ", "), "\n") diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index ce4877d2..0b647b81 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -1,22 +1,43 @@ -# Copyright (C) 2010-2019 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/> if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value_mean")) +#' Calculate the minimum error to assume in order to pass the variance test +#' +#' This function finds the smallest relative error still resulting in passing +#' the chi-squared test as defined in the FOCUS kinetics report from 2006. +#' +#' This function is used internally by \code{\link{summary.mkinfit}}. +#' +#' @param fit an object of class \code{\link{mkinfit}}. +#' @param alpha The confidence level chosen for the chi-squared test. +#' @importFrom stats qchisq aggregate +#' @return A dataframe with the following components: \item{err.min}{The +#' relative error, expressed as a fraction.} \item{n.optim}{The number of +#' optimised parameters attributed to the data series.} \item{df}{The number of +#' remaining degrees of freedom for the chi2 error level calculations. Note +#' that mean values are used for the chi2 statistic and therefore every time +#' point with observed values in the series only counts one time.} The +#' dataframe has one row for the total dataset and one further row for each +#' observed state variable in the model. +#' @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} +#' @keywords manip +#' @examples +#' +#' SFO_SFO = mkinmod(parent = mkinsub("SFO", to = "m1"), +#' m1 = mkinsub("SFO"), +#' use_of_ff = "max") +#' +#' fit_FOCUS_D = mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) +#' round(mkinerrmin(fit_FOCUS_D), 4) +#' \dontrun{ +#' fit_FOCUS_E = mkinfit(SFO_SFO, FOCUS_2006_E, quiet = TRUE) +#' round(mkinerrmin(fit_FOCUS_E), 4) +#' } +#' +#' @export mkinerrmin <- function(fit, alpha = 0.05) { parms.optim <- fit$par diff --git a/R/mkinerrplot.R b/R/mkinerrplot.R index 6153a3c0..36e22a43 100644 --- a/R/mkinerrplot.R +++ b/R/mkinerrplot.R @@ -1,22 +1,44 @@ -# Copyright (C) 2008-2014,2019 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/> if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual")) +#' Function to plot squared residuals and the error model for an mkin object +#' +#' This function plots the squared residuals for the specified subset of the +#' observed variables from an mkinfit object. In addition, one or more dashed +#' line(s) show the fitted error model. A combined plot of the fitted model +#' and this error model plot can be obtained with \code{\link{plot.mkinfit}} +#' using the argument \code{show_errplot = TRUE}. +#' +#' @param object A fit represented in an \code{\link{mkinfit}} object. +#' @param obs_vars A character vector of names of the observed variables for +#' which residuals should be plotted. Defaults to all observed variables in +#' the model +#' @param xlim plot range in x direction. +#' @param xlab Label for the x axis. +#' @param ylab Label for the y axis. +#' @param maxy Maximum value of the residuals. This is used for the scaling of +#' the y axis and defaults to "auto". +#' @param legend Should a legend be plotted? +#' @param lpos Where should the legend be placed? Default is "topright". Will +#' be passed on to \code{\link{legend}}. +#' @param col_obs Colors for the observed variables. +#' @param pch_obs Symbols to be used for the observed variables. +#' @param frame Should a frame be drawn around the plots? +#' @param \dots further arguments passed to \code{\link{plot}}. +#' @return Nothing is returned by this function, as it is called for its side +#' effect, namely to produce a plot. +#' @author Johannes Ranke +#' @seealso \code{\link{mkinplot}}, for a way to plot the data and the fitted +#' lines of the mkinfit object. +#' @keywords hplot +#' @examples +#' +#' \dontrun{ +#' model <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) +#' fit <- mkinfit(model, FOCUS_2006_D, error_model = "tc", quiet = TRUE) +#' mkinerrplot(fit) +#' } +#' +#' @export mkinerrplot <- function (object, obs_vars = names(object$mkinmod$map), xlim = c(0, 1.1 * max(object$data$predicted)), diff --git a/R/mkinfit.R b/R/mkinfit.R index d182f5d0..17fd59d0 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -1,904 +1,897 @@ -# Copyright (C) 2010-2019 Johannes Ranke
-# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG
-# 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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
-
-mkinfit <- function(mkinmod, observed,
- parms.ini = "auto",
- state.ini = "auto",
- err.ini = "auto",
- fixed_parms = NULL,
- fixed_initials = names(mkinmod$diffs)[-1],
- from_max_mean = FALSE,
- solution_type = c("auto", "analytical", "eigen", "deSolve"),
- method.ode = "lsoda",
- use_compiled = "auto",
- control = list(eval.max = 300, iter.max = 200),
- transform_rates = TRUE,
- transform_fractions = TRUE,
- quiet = FALSE,
- atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
- error_model = c("const", "obs", "tc"),
- error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", "OLS"),
- reweight.tol = 1e-8, reweight.max.iter = 10,
- trace_parms = FALSE,
- ...)
-{
- # Check mkinmod and generate a model for the variable whith the highest value
- # if a suitable string is given
- parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")
- if (class(mkinmod) != "mkinmod") {
- presumed_parent_name = observed[which.max(observed$value), "name"]
- if (mkinmod[[1]] %in% parent_models_available) {
- speclist <- list(list(type = mkinmod, sink = TRUE))
- names(speclist) <- presumed_parent_name
- mkinmod <- mkinmod(speclist = speclist)
- } else {
- stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ",
- paste(parent_models_available, collapse = ", "))
- }
- }
-
- # Get the names of the state variables in the model
- mod_vars <- names(mkinmod$diffs)
-
- # Get the names of observed variables
- obs_vars <- names(mkinmod$spec)
-
- # Subset observed data with names of observed data in the model and remove NA values
- observed <- subset(observed, name %in% obs_vars)
- observed <- subset(observed, !is.na(value))
-
- # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model)
- if (any(observed$value == 0)) {
- warning("Observations with value of zero were removed from the data")
- observed <- subset(observed, value != 0)
- }
-
- # Obtain data for decline from maximum mean value if requested
- if (from_max_mean) {
- # This is only used for simple decline models
- if (length(obs_vars) > 1)
- stop("Decline from maximum is only implemented for models with a single observed variable")
-
- means <- aggregate(value ~ time, data = observed, mean, na.rm=TRUE)
- t_of_max <- means[which.max(means$value), "time"]
- observed <- subset(observed, time >= t_of_max)
- observed$time <- observed$time - t_of_max
- }
-
- # Number observations used for fitting
- n_observed <- nrow(observed)
-
- # Define starting values for parameters where not specified by the user
- if (parms.ini[[1]] == "auto") parms.ini = vector()
-
- # Warn for inital parameter specifications that are not in the model
- wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms)
- if (length(wrongpar.names) > 0) {
- warning("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "),
- " not used in the model")
- parms.ini <- parms.ini[setdiff(names(parms.ini), wrongpar.names)]
- }
-
- # Warn that the sum of formation fractions may exceed one if they are not
- # fitted in the transformed way
- if (mkinmod$use_of_ff == "max" & transform_fractions == FALSE) {
- warning("The sum of formation fractions may exceed one if you do not use ",
- "transform_fractions = TRUE." )
- for (box in mod_vars) {
- # Stop if formation fractions are not transformed and we have no sink
- if (mkinmod$spec[[box]]$sink == FALSE) {
- stop("If formation fractions are not transformed during the fitting, ",
- "it is not supported to turn off pathways to sink.\n ",
- "Consider turning on the transformation of formation fractions or ",
- "setting up a model with use_of_ff = 'min'.\n")
- }
- }
- }
-
- # Do not allow fixing formation fractions if we are using the ilr transformation,
- # this is not supported
- if (transform_fractions == TRUE && length(fixed_parms) > 0) {
- if (any(grepl("^f_", fixed_parms))) {
- stop("Fixing formation fractions is not supported when using the ilr ",
- "transformation.")
- }
- }
-
- # Set initial parameter values, including a small increment (salt)
- # to avoid linear dependencies (singular matrix) in Eigenvalue based solutions
- k_salt = 0
- defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
- for (parmname in defaultpar.names) {
- # Default values for rate constants, depending on the parameterisation
- if (grepl("^k", parmname)) {
- parms.ini[parmname] = 0.1 + k_salt
- k_salt = k_salt + 1e-4
- }
- # Default values for rate constants for reversible binding
- if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1
- if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02
- # Default values for IORE exponents
- if (grepl("^N", parmname)) parms.ini[parmname] = 1.1
- # Default values for the FOMC, DFOP and HS models
- if (parmname == "alpha") parms.ini[parmname] = 1
- if (parmname == "beta") parms.ini[parmname] = 10
- if (parmname == "k1") parms.ini[parmname] = 0.1
- if (parmname == "k2") parms.ini[parmname] = 0.01
- if (parmname == "tb") parms.ini[parmname] = 5
- if (parmname == "g") parms.ini[parmname] = 0.5
- if (parmname == "kmax") parms.ini[parmname] = 0.1
- if (parmname == "k0") parms.ini[parmname] = 0.0001
- if (parmname == "r") parms.ini[parmname] = 0.2
- }
- # Default values for formation fractions in case they are present
- for (box in mod_vars) {
- f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)]
- if (length(f_names) > 0) {
- # We need to differentiate between default and specified fractions
- # and set the unspecified to 1 - sum(specified)/n_unspecified
- f_default_names <- intersect(f_names, defaultpar.names)
- f_specified_names <- setdiff(f_names, defaultpar.names)
- sum_f_specified = sum(parms.ini[f_specified_names])
- if (sum_f_specified > 1) {
- stop("Starting values for the formation fractions originating from ",
- box, " sum up to more than 1.")
- }
- if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1
- else {
- n_unspecified = length(f_default_names)
- }
- parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified
- }
- }
-
- # Set default for state.ini if appropriate
- parent_name = names(mkinmod$spec)[[1]]
- if (state.ini[1] == "auto") {
- parent_time_0 = subset(observed, time == 0 & name == parent_name)$value
- parent_time_0_mean = mean(parent_time_0, na.rm = TRUE)
- if (is.na(parent_time_0_mean)) {
- state.ini = c(100, rep(0, length(mkinmod$diffs) - 1))
- } else {
- state.ini = c(parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1))
- }
- }
-
- # Name the inital state variable values if they are not named yet
- if(is.null(names(state.ini))) names(state.ini) <- mod_vars
-
- # Transform initial parameter values for fitting
- transparms.ini <- transform_odeparms(parms.ini, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Parameters to be optimised:
- # Kinetic parameters in parms.ini whose names are not in fixed_parms
- parms.fixed <- parms.ini[fixed_parms]
- parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)]
-
- transparms.fixed <- transform_odeparms(parms.fixed, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- transparms.optim <- transform_odeparms(parms.optim, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Inital state variables in state.ini whose names are not in fixed_initials
- state.ini.fixed <- state.ini[fixed_initials]
- state.ini.optim <- state.ini[setdiff(names(state.ini), fixed_initials)]
-
- # Preserve names of state variables before renaming initial state variable
- # parameters
- state.ini.optim.boxnames <- names(state.ini.optim)
- state.ini.fixed.boxnames <- names(state.ini.fixed)
- if(length(state.ini.optim) > 0) {
- names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
- }
- if(length(state.ini.fixed) > 0) {
- names(state.ini.fixed) <- paste(names(state.ini.fixed), "0", sep="_")
- }
-
- # Decide if the solution of the model can be based on a simple analytical
- # formula, the spectral decomposition of the matrix (fundamental system)
- # or a numeric ode solver from the deSolve package
- # Prefer deSolve over eigen if a compiled model is present and use_compiled
- # is not set to FALSE
- solution_type = match.arg(solution_type)
- if (solution_type == "analytical" && length(mkinmod$spec) > 1)
- stop("Analytical solution not implemented for models with metabolites.")
- if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
- stop("Eigenvalue based solution not possible, coefficient matrix not present.")
- if (solution_type == "auto") {
- if (length(mkinmod$spec) == 1) {
- solution_type = "analytical"
- } else {
- if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) {
- solution_type = "deSolve"
- } else {
- if (is.matrix(mkinmod$coefmat)) {
- solution_type = "eigen"
- if (max(observed$value, na.rm = TRUE) < 0.1) {
- stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone")
- }
- } else {
- solution_type = "deSolve"
- }
- }
- }
- }
-
- # Get the error model and the algorithm for fitting
- err_mod <- match.arg(error_model)
- error_model_algorithm = match.arg(error_model_algorithm)
- if (error_model_algorithm == "OLS") {
- if (err_mod != "const") stop("OLS is only appropriate for constant variance")
- }
- if (error_model_algorithm == "auto") {
- error_model_algorithm = switch(err_mod,
- const = "OLS", obs = "d_3", tc = "d_3")
- }
- errparm_names <- switch(err_mod,
- "const" = "sigma",
- "obs" = paste0("sigma_", obs_vars),
- "tc" = c("sigma_low", "rsd_high"))
- errparm_names_optim <- if (error_model_algorithm == "OLS") NULL else errparm_names
-
- # Define starting values for the error model
- if (err.ini[1] != "auto") {
- if (!identical(names(err.ini), errparm_names)) {
- stop("Please supply initial values for error model components ", paste(errparm_names, collapse = ", "))
- } else {
- errparms = err.ini
- }
- } else {
- if (err_mod == "const") {
- errparms = 3
- }
- if (err_mod == "obs") {
- errparms = rep(3, length(obs_vars))
- }
- if (err_mod == "tc") {
- errparms <- c(sigma_low = 0.1, rsd_high = 0.1)
- }
- names(errparms) <- errparm_names
- }
- if (error_model_algorithm == "OLS") {
- errparms_optim <- NULL
- } else {
- errparms_optim <- errparms
- }
-
- # Define outtimes for model solution.
- # Include time points at which observed data are available
- outtimes = sort(unique(c(observed$time, seq(min(observed$time),
- max(observed$time),
- length.out = n.outtimes))))
-
- # Define the objective function for optimisation, including (back)transformations
- cost_function <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...)
- {
- assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter
-
- # Trace parameter values if requested and if we are actually optimising
- if(trace_parms & update_data) cat(P, "\n")
-
- # Determine local parameter values for the cost estimation
- if (is.numeric(fixed_degparms)) {
- cost_degparms <- fixed_degparms
- cost_errparms <- P
- degparms_fixed = TRUE
- } else {
- degparms_fixed = FALSE
- }
-
- if (is.numeric(fixed_errparms)) {
- cost_degparms <- P
- cost_errparms <- fixed_errparms
- errparms_fixed = TRUE
- } else {
- errparms_fixed = FALSE
- }
-
- if (OLS) {
- cost_degparms <- P
- cost_errparms <- numeric(0)
- }
-
- if (!OLS & !degparms_fixed & !errparms_fixed) {
- cost_degparms <- P[1:(length(P) - length(errparms))]
- cost_errparms <- P[(length(cost_degparms) + 1):length(P)]
- }
-
- # Initial states for t0
- if(length(state.ini.optim) > 0) {
- odeini <- c(cost_degparms[1:length(state.ini.optim)], state.ini.fixed)
- names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
- } else {
- odeini <- state.ini.fixed
- names(odeini) <- state.ini.fixed.boxnames
- }
-
- odeparms.optim <- cost_degparms[(length(state.ini.optim) + 1):length(cost_degparms)]
-
- if (trans == TRUE) {
- odeparms <- c(odeparms.optim, transparms.fixed)
- parms <- backtransform_odeparms(odeparms, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- } else {
- parms <- c(odeparms.optim, parms.fixed)
- }
-
- # Solve the system with current parameter values
- out <- mkinpredict(mkinmod, parms,
- odeini, outtimes,
- solution_type = solution_type,
- use_compiled = use_compiled,
- method.ode = method.ode,
- atol = atol, rtol = rtol, ...)
-
- out_long <- mkin_wide_to_long(out, time = "time")
-
- if (err_mod == "const") {
- observed$std <- if (OLS) NA else cost_errparms["sigma"]
- }
- if (err_mod == "obs") {
- std_names <- paste0("sigma_", observed$name)
- observed$std <- cost_errparms[std_names]
- }
- if (err_mod == "tc") {
- tmp <- merge(observed, out_long, by = c("time", "name"))
- tmp$name <- ordered(tmp$name, levels = obs_vars)
- tmp <- tmp[order(tmp$name, tmp$time), ]
- observed$std <- sqrt(cost_errparms["sigma_low"]^2 + tmp$value.y^2 * cost_errparms["rsd_high"]^2)
- }
-
- cost_data <- merge(observed[c("name", "time", "value", "std")], out_long,
- by = c("name", "time"), suffixes = c(".observed", ".predicted"))
-
- if (OLS) {
- # Cost is the sum of squared residuals
- cost <- with(cost_data, sum((value.observed - value.predicted)^2))
- } else {
- # Cost is the negative log-likelihood
- cost <- - with(cost_data,
- sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE)))
- }
-
- # We update the current cost and data during the optimisation, not
- # during hessian calculations
- if (update_data) {
-
- assign("out_predicted", out_long, inherits = TRUE)
- assign("current_data", cost_data, inherits = TRUE)
-
- if (cost < cost.current) {
- assign("cost.current", cost, inherits = TRUE)
- if (!quiet) cat(ifelse(OLS, "Sum of squared residuals", "Negative log-likelihood"),
- " at call ", calls, ": ", cost.current, "\n", sep = "")
- }
- }
- return(cost)
- }
-
- names_optim <- c(names(state.ini.optim),
- names(transparms.optim),
- errparm_names_optim)
- n_optim <- length(names_optim)
-
- # Define lower and upper bounds other than -Inf and Inf for parameters
- # for which no internal transformation is requested in the call to mkinfit
- # and for optimised error model parameters
- lower <- rep(-Inf, n_optim)
- upper <- rep(Inf, n_optim)
- names(lower) <- names(upper) <- names_optim
-
- # IORE exponents are not transformed, but need a lower bound
- index_N <- grep("^N", names(lower))
- lower[index_N] <- 0
-
- if (!transform_rates) {
- index_k <- grep("^k_", names(lower))
- lower[index_k] <- 0
- index_k__iore <- grep("^k__iore_", names(lower))
- lower[index_k__iore] <- 0
- other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb", "r"), names(lower))
- lower[other_rate_parms] <- 0
- }
-
- if (!transform_fractions) {
- index_f <- grep("^f_", names(upper))
- lower[index_f] <- 0
- upper[index_f] <- 1
- other_fraction_parms <- intersect(c("g"), names(upper))
- lower[other_fraction_parms] <- 0
- upper[other_fraction_parms] <- 1
- }
-
- if (err_mod == "const") {
- if (error_model_algorithm != "OLS") {
- lower["sigma"] <- 0
- }
- }
- if (err_mod == "obs") {
- index_sigma <- grep("^sigma_", names(lower))
- lower[index_sigma] <- 0
- }
- if (err_mod == "tc") {
- lower["sigma_low"] <- 0
- lower["rsd_high"] <- 0
- }
-
- # Counter for cost function evaluations
- calls = 0
- cost.current <- Inf
- out_predicted <- NA
- current_data <- NA
-
- # Show parameter names if tracing is requested
- if(trace_parms) cat(names_optim, "\n")
-
- # browser()
-
- # Do the fit and take the time until the hessians are calculated
- fit_time <- system.time({
- degparms <- c(state.ini.optim, transparms.optim)
- n_degparms <- length(degparms)
- degparms_index <- seq(1, n_degparms)
- errparms_index <- seq(n_degparms + 1, length.out = length(errparms))
-
- if (error_model_algorithm == "d_3") {
- if (!quiet) message("Directly optimising the complete model")
- parms.start <- c(degparms, errparms)
- fit_direct <- nlminb(parms.start, cost_function,
- lower = lower[names(parms.start)],
- upper = upper[names(parms.start)],
- control = control, ...)
- fit_direct$logLik <- - cost.current
- if (error_model_algorithm == "direct") {
- degparms <- fit_direct$par[degparms_index]
- errparms <- fit_direct$par[errparms_index]
- } else {
- cost.current <- Inf # reset to avoid conflict with the OLS step
- }
- }
- if (error_model_algorithm != "direct") {
- if (!quiet) message("Ordinary least squares optimisation")
- fit <- nlminb(degparms, cost_function, control = control,
- lower = lower[names(degparms)],
- upper = upper[names(degparms)], OLS = TRUE, ...)
- degparms <- fit$par
-
- # Get the maximum likelihood estimate for sigma at the optimum parameter values
- current_data$residual <- current_data$value.observed - current_data$value.predicted
- sigma_mle <- sqrt(sum(current_data$residual^2)/nrow(current_data))
-
- # Use that estimate for the constant variance, or as first guess if err_mod = "obs"
- if (err_mod != "tc") {
- errparms[names(errparms)] <- sigma_mle
- }
- fit$par <- c(fit$par, errparms)
-
- cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
- fit$logLik <- - cost.current
- }
- if (error_model_algorithm %in% c("threestep", "fourstep", "d_3")) {
- if (!quiet) message("Optimising the error model")
- fit <- nlminb(errparms, cost_function, control = control,
- lower = lower[names(errparms)],
- upper = upper[names(errparms)],
- fixed_degparms = degparms, ...)
- errparms <- fit$par
- }
- if (error_model_algorithm == "fourstep") {
- if (!quiet) message("Optimising the degradation model")
- fit <- nlminb(degparms, cost_function, control = control,
- lower = lower[names(degparms)],
- upper = upper[names(degparms)],
- fixed_errparms = errparms, ...)
- degparms <- fit$par
- }
- if (error_model_algorithm %in%
- c("direct", "twostep", "threestep", "fourstep", "d_3")) {
- if (!quiet) message("Optimising the complete model")
- parms.start <- c(degparms, errparms)
- fit <- nlminb(parms.start, cost_function,
- lower = lower[names(parms.start)],
- upper = upper[names(parms.start)],
- control = control, ...)
- degparms <- fit$par[degparms_index]
- errparms <- fit$par[errparms_index]
- fit$logLik <- - cost.current
-
- if (error_model_algorithm == "d_3") {
- d_3_messages = c(
- same = "Direct fitting and three-step fitting yield approximately the same likelihood",
- threestep = "Three-step fitting yielded a higher likelihood than direct fitting",
- direct = "Direct fitting yielded a higher likelihood than three-step fitting")
- rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
- if (rel_diff < 0.0001) {
- if (!quiet) message(d_3_messages["same"])
- fit$d_3_message <- d_3_messages["same"]
- } else {
- if (fit$logLik > fit_direct$logLik) {
- if (!quiet) message(d_3_messages["threestep"])
- fit$d_3_message <- d_3_messages["threestep"]
- } else {
- if (!quiet) message(d_3_messages["direct"])
- fit <- fit_direct
- fit$d_3_message <- d_3_messages["direct"]
- }
- }
- }
- }
- if (err_mod != "const" & error_model_algorithm == "IRLS") {
- reweight.diff <- 1
- n.iter <- 0
- errparms_last <- errparms
-
- while (reweight.diff > reweight.tol &
- n.iter < reweight.max.iter) {
-
- if (!quiet) message("Optimising the error model")
- fit <- nlminb(errparms, cost_function, control = control,
- lower = lower[names(errparms)],
- upper = upper[names(errparms)],
- fixed_degparms = degparms, ...)
- errparms <- fit$par
-
- if (!quiet) message("Optimising the degradation model")
- fit <- nlminb(degparms, cost_function, control = control,
- lower = lower[names(degparms)],
- upper = upper[names(degparms)],
- fixed_errparms = errparms, ...)
- degparms <- fit$par
-
- reweight.diff <- dist(rbind(errparms, errparms_last))
- errparms_last <- errparms
-
- fit$par <- c(fit$par, errparms)
- cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
- fit$logLik <- - cost.current
- }
- }
-
- fit$hessian <- try(numDeriv::hessian(cost_function, c(degparms, errparms), OLS = FALSE,
- update_data = FALSE), silent = TRUE)
-
- # Backtransform parameters
- bparms.optim = backtransform_odeparms(fit$par, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- bparms.fixed = c(state.ini.fixed, parms.fixed)
- bparms.all = c(bparms.optim, parms.fixed)
-
- fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.all, errparms),
- OLS = FALSE, trans = FALSE, update_data = FALSE), silent = TRUE)
- })
-
- fit$error_model_algorithm <- error_model_algorithm
-
- if (fit$convergence != 0) {
- fit$warning = paste0("Optimisation did not converge:\n", fit$message)
- warning(fit$warning)
- } else {
- if(!quiet) message("Optimisation successfully terminated.\n")
- }
-
- # We need to return some more data for summary and plotting
- fit$solution_type <- solution_type
- fit$transform_rates <- transform_rates
- fit$transform_fractions <- transform_fractions
- fit$reweight.tol <- reweight.tol
- fit$reweight.max.iter <- reweight.max.iter
- fit$control <- control
- fit$calls <- calls
- fit$time <- fit_time
-
- # We also need the model for summary and plotting
- fit$mkinmod <- mkinmod
-
- # We need data and predictions for summary and plotting
- fit$observed <- observed
- fit$obs_vars <- obs_vars
- fit$predicted <- out_predicted
-
- # Residual sum of squares as a function of the fitted parameters
- fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE)
-
- # Log-likelihood with possibility to fix degparms or errparms
- fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE) {
- - cost_function(P, fixed_degparms = fixed_degparms,
- fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE)
- }
-
- # Collect initial parameter values in three dataframes
- fit$start <- data.frame(value = c(state.ini.optim,
- parms.optim, errparms_optim))
- fit$start$type = c(rep("state", length(state.ini.optim)),
- rep("deparm", length(parms.optim)),
- rep("error", length(errparms_optim)))
-
- fit$start_transformed = data.frame(
- value = c(state.ini.optim, transparms.optim, errparms_optim),
- lower = lower,
- upper = upper)
-
- fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))
- fit$fixed$type = c(rep("state", length(state.ini.fixed)),
- rep("deparm", length(parms.fixed)))
-
- # Sort observed, predicted and residuals
- current_data$name <- ordered(current_data$name, levels = obs_vars)
-
- ordered_data <- current_data[order(current_data$name, current_data$time), ]
-
- fit$data <- data.frame(time = ordered_data$time,
- variable = ordered_data$name,
- observed = ordered_data$value.observed,
- predicted = ordered_data$value.predicted)
-
- fit$data$residual <- fit$data$observed - fit$data$predicted
-
- fit$atol <- atol
- fit$rtol <- rtol
- fit$err_mod <- err_mod
-
- # Return different sets of backtransformed parameters for summary and plotting
- fit$bparms.optim <- bparms.optim
- fit$bparms.fixed <- bparms.fixed
-
- # Return ode and state parameters for further fitting
- fit$bparms.ode <- bparms.all[mkinmod$parms]
- fit$bparms.state <- c(bparms.all[setdiff(names(bparms.all), names(fit$bparms.ode))],
- state.ini.fixed)
- names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state))
-
- fit$errparms <- errparms
- fit$df.residual <- n_observed - length(c(degparms, errparms))
-
- fit$date <- date()
- fit$version <- as.character(utils::packageVersion("mkin"))
- fit$Rversion <- paste(R.version$major, R.version$minor, sep=".")
-
- class(fit) <- c("mkinfit", "modFit")
- return(fit)
-}
-
-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
- param <- object$par
- pnames <- names(param)
- bpnames <- names(object$bparms.optim)
- epnames <- names(object$errparms)
- p <- length(param)
- mod_vars <- names(object$mkinmod$diffs)
- covar <- try(solve(object$hessian), silent = TRUE)
- covar_notrans <- try(solve(object$hessian_notrans), silent = TRUE)
- rdf <- object$df.residual
-
- if (!is.numeric(covar) | is.na(covar[1])) {
- covar <- NULL
- se <- lci <- uci <- rep(NA, p)
- } else {
- rownames(covar) <- colnames(covar) <- pnames
- se <- sqrt(diag(covar))
- lci <- param + qt(alpha/2, rdf) * se
- uci <- param + qt(1-alpha/2, rdf) * se
- }
-
- beparms.optim <- c(object$bparms.optim, object$par[epnames])
- if (!is.numeric(covar_notrans) | is.na(covar_notrans[1])) {
- covar_notrans <- NULL
- se_notrans <- tval <- pval <- rep(NA, p)
- } else {
- rownames(covar_notrans) <- colnames(covar_notrans) <- c(bpnames, epnames)
- se_notrans <- sqrt(diag(covar_notrans))
- tval <- beparms.optim / se_notrans
- pval <- pt(abs(tval), rdf, lower.tail = FALSE)
- }
-
- names(se) <- pnames
-
- param <- cbind(param, se, lci, uci)
- dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))
-
- bparam <- cbind(Estimate = beparms.optim, se_notrans,
- "t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA)
-
- # Transform boundaries of CI for one parameter at a time,
- # with the exception of sets of formation fractions (single fractions are OK).
- f_names_skip <- character(0)
- for (box in mod_vars) { # Figure out sets of fractions to skip
- f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
- n_paths <- length(f_names)
- if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
- }
-
- for (pname in pnames) {
- if (!pname %in% f_names_skip) {
- par.lower <- param[pname, "Lower"]
- par.upper <- param[pname, "Upper"]
- names(par.lower) <- names(par.upper) <- pname
- bpl <- backtransform_odeparms(par.lower, object$mkinmod,
- object$transform_rates,
- object$transform_fractions)
- bpu <- backtransform_odeparms(par.upper, object$mkinmod,
- object$transform_rates,
- object$transform_fractions)
- bparam[names(bpl), "Lower"] <- bpl
- bparam[names(bpu), "Upper"] <- bpu
- }
- }
- bparam[epnames, c("Lower", "Upper")] <- param[epnames, c("Lower", "Upper")]
-
- ans <- list(
- version = as.character(utils::packageVersion("mkin")),
- Rversion = paste(R.version$major, R.version$minor, sep="."),
- date.fit = object$date,
- date.summary = date(),
- solution_type = object$solution_type,
- warning = object$warning,
- use_of_ff = object$mkinmod$use_of_ff,
- error_model_algorithm = object$error_model_algorithm,
- df = c(p, rdf),
- covar = covar,
- covar_notrans = covar_notrans,
- err_mod = object$err_mod,
- niter = object$iterations,
- calls = object$calls,
- time = object$time,
- par = param,
- bpar = bparam)
-
- if (!is.null(object$version)) {
- ans$fit_version <- object$version
- ans$fit_Rversion <- object$Rversion
- }
-
- ans$diffs <- object$mkinmod$diffs
- if(data) ans$data <- object$data
- ans$start <- object$start
- ans$start_transformed <- object$start_transformed
-
- ans$fixed <- object$fixed
-
- ans$errmin <- mkinerrmin(object, alpha = 0.05)
-
- if (object$calls > 0) {
- if (!is.null(ans$covar)){
- Corr <- cov2cor(ans$covar)
- rownames(Corr) <- colnames(Corr) <- rownames(ans$par)
- ans$Corr <- Corr
- } else {
- warning("Could not calculate correlation; no covariance matrix")
- }
- }
-
- ans$bparms.ode <- object$bparms.ode
- ep <- endpoints(object)
- if (length(ep$ff) != 0)
- ans$ff <- ep$ff
- if (distimes) ans$distimes <- ep$distimes
- if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
- if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message
- class(ans) <- c("summary.mkinfit", "summary.modFit")
- return(ans)
-}
-
-# Expanded from print.summary.modFit
-print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
- if (is.null(x$fit_version)) {
- cat("mkin version: ", x$version, "\n")
- cat("R version: ", x$Rversion, "\n")
- } else {
- cat("mkin version used for fitting: ", x$fit_version, "\n")
- cat("R version used for fitting: ", x$fit_Rversion, "\n")
- }
-
- cat("Date of fit: ", x$date.fit, "\n")
- cat("Date of summary:", x$date.summary, "\n")
-
- if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n")
-
- cat("\nEquations:\n")
- nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
- writeLines(strwrap(nice_diffs, exdent = 11))
- df <- x$df
- rdf <- df[2]
-
- cat("\nModel predictions using solution type", x$solution_type, "\n")
-
- cat("\nFitted using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n")
-
- if (!is.null(x$err_mod)) {
- cat("\nError model: ")
- cat(switch(x$err_mod,
- const = "Constant variance",
- obs = "Variance unique to each observed variable",
- tc = "Two-component variance function"), "\n")
-
- cat("\nError model algorithm:", x$error_model_algorithm, "\n")
- if (!is.null(x$d_3_message)) cat(x$d_3_message, "\n")
- }
-
- cat("\nStarting values for parameters to be optimised:\n")
- print(x$start)
-
- cat("\nStarting values for the transformed parameters actually optimised:\n")
- print(x$start_transformed)
-
- cat("\nFixed parameter values:\n")
- if(length(x$fixed$value) == 0) cat("None\n")
- else print(x$fixed)
-
- cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
- print(signif(x$par, digits = digits))
-
- if (x$calls > 0) {
- cat("\nParameter correlation:\n")
- if (!is.null(x$covar)){
- print(x$Corr, digits = digits, ...)
- } else {
- cat("No covariance matrix")
- }
- }
-
- cat("\nBacktransformed parameters:\n")
- cat("Confidence intervals for internally transformed parameters are asymmetric.\n")
- if ((x$version) < "0.9-36") {
- cat("To get the usual (questionable) t-test, upgrade mkin and repeat the fit.\n")
- print(signif(x$bpar, digits = digits))
- } else {
- cat("t-test (unrealistically) based on the assumption of normal distribution\n")
- cat("for estimators of untransformed parameters.\n")
- print(signif(x$bpar[, c(1, 3, 4, 5, 6)], digits = digits))
- }
-
- cat("\nFOCUS Chi2 error levels in percent:\n")
- x$errmin$err.min <- 100 * x$errmin$err.min
- print(x$errmin, digits=digits,...)
-
- printSFORB <- !is.null(x$SFORB)
- if(printSFORB){
- cat("\nEstimated Eigenvalues of SFORB model(s):\n")
- print(x$SFORB, digits=digits,...)
- }
-
- printff <- !is.null(x$ff)
- if(printff){
- cat("\nResulting formation fractions:\n")
- print(data.frame(ff = x$ff), digits=digits,...)
- }
-
- printdistimes <- !is.null(x$distimes)
- if(printdistimes){
- cat("\nEstimated disappearance times:\n")
- print(x$distimes, digits=digits,...)
- }
-
- printdata <- !is.null(x$data)
- if (printdata){
- cat("\nData:\n")
- print(format(x$data, digits = digits, ...), row.names = FALSE)
- }
-
- invisible(x)
-}
-# vim: set ts=2 sw=2 expandtab:
+if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) + +#' Fit a kinetic model to data with one or more state variables +#' +#' This function maximises the likelihood of the observed data using the Port +#' algorithm \code{\link{nlminb}}, and the specified initial or fixed +#' parameters and starting values. In each step of the optimsation, the +#' kinetic model is solved using the function \code{\link{mkinpredict}}. The +#' parameters of the selected error model are fitted simultaneously with the +#' degradation model parameters, as both of them are arguments of the +#' likelihood function. +#' +#' Per default, parameters in the kinetic models are internally transformed in +#' order to better satisfy the assumption of a normal distribution of their +#' estimators. +#' +#' @param mkinmod A list of class \code{\link{mkinmod}}, containing the kinetic +#' model to be fitted to the data, or one of the shorthand names ("SFO", +#' "FOMC", "DFOP", "HS", "SFORB", "IORE"). If a shorthand name is given, a +#' parent only degradation model is generated for the variable with the +#' highest value in \code{observed}. +#' @param observed A dataframe with the observed data. The first column called +#' "name" must contain the name of the observed variable for each data point. +#' The second column must contain the times of observation, named "time". +#' The third column must be named "value" and contain the observed values. +#' Zero values in the "value" column will be removed, with a warning, in +#' order to avoid problems with fitting the two-component error model. This +#' is not expected to be a problem, because in general, values of zero are +#' not observed in degradation data, because there is a lower limit of +#' detection. +#' @param parms.ini A named vector of initial values for the parameters, +#' including parameters to be optimised and potentially also fixed parameters +#' as indicated by \code{fixed_parms}. If set to "auto", initial values for +#' rate constants are set to default values. Using parameter names that are +#' not in the model gives an error. +#' +#' It is possible to only specify a subset of the parameters that the model +#' needs. You can use the parameter lists "bparms.ode" from a previously +#' fitted model, which contains the differential equation parameters from +#' this model. This works nicely if the models are nested. An example is +#' given below. +#' @param state.ini A named vector of initial values for the state variables of +#' the model. In case the observed variables are represented by more than one +#' model variable, the names will differ from the names of the observed +#' variables (see \code{map} component of \code{\link{mkinmod}}). The default +#' is to set the initial value of the first model variable to the mean of the +#' time zero values for the variable with the maximum observed value, and all +#' others to 0. If this variable has no time zero observations, its initial +#' value is set to 100. +#' @param err.ini A named vector of initial values for the error model +#' parameters to be optimised. If set to "auto", initial values are set to +#' default values. Otherwise, inital values for all error model parameters +#' must be given. +#' @param fixed_parms The names of parameters that should not be optimised but +#' rather kept at the values specified in \code{parms.ini}. +#' @param fixed_initials The names of model variables for which the initial +#' state at time 0 should be excluded from the optimisation. Defaults to all +#' state variables except for the first one. +#' @param from_max_mean If this is set to TRUE, and the model has only one +#' observed variable, then data before the time of the maximum observed value +#' (after averaging for each sampling time) are discarded, and this time is +#' subtracted from all remaining time values, so the time of the maximum +#' observed mean value is the new time zero. +#' @param solution_type If set to "eigen", the solution of the system of +#' differential equations is based on the spectral decomposition of the +#' coefficient matrix in cases that this is possible. If set to "deSolve", a +#' numerical ode solver from package \code{\link{deSolve}} is used. If set to +#' "analytical", an analytical solution of the model is used. This is only +#' implemented for simple degradation experiments with only one state +#' variable, i.e. with no metabolites. The default is "auto", which uses +#' "analytical" if possible, otherwise "deSolve" if a compiler is present, +#' and "eigen" if no compiler is present and the model can be expressed using +#' eigenvalues and eigenvectors. This argument is passed on to the helper +#' function \code{\link{mkinpredict}}. +#' @param method.ode The solution method passed via \code{\link{mkinpredict}} +#' to \code{\link{ode}} in case the solution type is "deSolve". The default +#' "lsoda" is performant, but sometimes fails to converge. +#' @param use_compiled If set to \code{FALSE}, no compiled version of the +#' \code{\link{mkinmod}} model is used in the calls to +#' \code{\link{mkinpredict}} even if a compiled version is present. +#' @param control A list of control arguments passed to \code{\link{nlminb}}. +#' @param transform_rates Boolean specifying if kinetic rate constants should +#' be transformed in the model specification used in the fitting for better +#' compliance with the assumption of normal distribution of the estimator. If +#' TRUE, also alpha and beta parameters of the FOMC model are +#' log-transformed, as well as k1 and k2 rate constants for the DFOP and HS +#' models and the break point tb of the HS model. If FALSE, zero is used as +#' a lower bound for the rates in the optimisation. +#' @param transform_fractions Boolean specifying if formation fractions +#' constants should be transformed in the model specification used in the +#' fitting for better compliance with the assumption of normal distribution +#' of the estimator. The default (TRUE) is to do transformations. If TRUE, +#' the g parameter of the DFOP and HS models are also transformed, as they +#' can also be seen as compositional data. The transformation used for these +#' transformations is the \code{\link{ilr}} transformation. +#' @param quiet Suppress printing out the current value of the negative +#' log-likelihood after each improvement? +#' @param atol Absolute error tolerance, passed to \code{\link{ode}}. Default +#' is 1e-8, lower than in \code{\link{lsoda}}. +#' @param rtol Absolute error tolerance, passed to \code{\link{ode}}. Default +#' is 1e-10, much lower than in \code{\link{lsoda}}. +#' @param n.outtimes The length of the dataseries that is produced by the model +#' prediction function \code{\link{mkinpredict}}. This impacts the accuracy +#' of the numerical solver if that is used (see \code{solution_type} +#' argument. The default value is 100. +#' @param error_model If the error model is "const", a constant standard +#' deviation is assumed. +#' +#' If the error model is "obs", each observed variable is assumed to have its +#' own variance. +#' +#' If the error model is "tc" (two-component error model), a two component +#' error model similar to the one described by Rocke and Lorenzato (1995) is +#' used for setting up the likelihood function. Note that this model +#' deviates from the model by Rocke and Lorenzato, as their model implies +#' that the errors follow a lognormal distribution for large values, not a +#' normal distribution as assumed by this method. +#' @param error_model_algorithm If "auto", the selected algorithm depends on +#' the error model. If the error model is "const", unweighted nonlinear +#' least squares fitting ("OLS") is selected. If the error model is "obs", or +#' "tc", the "d_3" algorithm is selected. +#' +#' The algorithm "d_3" will directly minimize the negative log-likelihood and +#' - independently - also use the three step algorithm described below. The +#' fit with the higher likelihood is returned. +#' +#' The algorithm "direct" will directly minimize the negative log-likelihood. +#' +#' The algorithm "twostep" will minimize the negative log-likelihood after an +#' initial unweighted least squares optimisation step. +#' +#' The algorithm "threestep" starts with unweighted least squares, then +#' optimizes only the error model using the degradation model parameters +#' found, and then minimizes the negative log-likelihood with free +#' degradation and error model parameters. +#' +#' The algorithm "fourstep" starts with unweighted least squares, then +#' optimizes only the error model using the degradation model parameters +#' found, then optimizes the degradation model again with fixed error model +#' parameters, and finally minimizes the negative log-likelihood with free +#' degradation and error model parameters. +#' +#' The algorithm "IRLS" (Iteratively Reweighted Least Squares) starts with +#' unweighted least squares, and then iterates optimization of the error +#' model parameters and subsequent optimization of the degradation model +#' using those error model parameters, until the error model parameters +#' converge. +#' @param reweight.tol Tolerance for the convergence criterion calculated from +#' the error model parameters in IRLS fits. +#' @param reweight.max.iter Maximum number of iterations in IRLS fits. +#' @param trace_parms Should a trace of the parameter values be listed? +#' @param \dots Further arguments that will be passed on to +#' \code{\link{deSolve}}. +#' @importFrom stats nlminb aggregate dist +#' @return A list with "mkinfit" in the class attribute. A summary can be +#' obtained by \code{\link{summary.mkinfit}}. +#' @note When using the "IORE" submodel for metabolites, fitting with +#' "transform_rates = TRUE" (the default) often leads to failures of the +#' numerical ODE solver. In this situation it may help to switch off the +#' internal rate transformation. +#' @author Johannes Ranke +#' @seealso Plotting methods \code{\link{plot.mkinfit}} and +#' \code{\link{mkinparplot}}. +#' +#' Comparisons of models fitted to the same data can be made using +#' \code{\link{AIC}} by virtue of the method \code{\link{logLik.mkinfit}}. +#' +#' Fitting of several models to several datasets in a single call to +#' \code{\link{mmkin}}. +#' @source Rocke, David M. und Lorenzato, Stefan (1995) A two-component model +#' for measurement error in analytical chemistry. Technometrics 37(2), 176-184. +#' @examples +#' +#' # Use shorthand notation for parent only degradation +#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) +#' summary(fit) +#' +#' # One parent compound, one metabolite, both single first order. +#' # Use mkinsub for convenience in model formulation. Pathway to sink included per default. +#' SFO_SFO <- mkinmod( +#' parent = mkinsub("SFO", "m1"), +#' m1 = mkinsub("SFO")) +#' # Fit the model to the FOCUS example dataset D using defaults +#' print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D, +#' solution_type = "eigen", quiet = TRUE))) +#' coef(fit) +#' endpoints(fit) +#' \dontrun{ +#' # deSolve is slower when no C compiler (gcc) was available during model generation +#' print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, +#' solution_type = "deSolve"))) +#' coef(fit.deSolve) +#' endpoints(fit.deSolve) +#' } +#' +#' # Use stepwise fitting, using optimised parameters from parent only fit, FOMC +#' \dontrun{ +#' FOMC_SFO <- mkinmod( +#' parent = mkinsub("FOMC", "m1"), +#' m1 = mkinsub("SFO")) +#' # Fit the model to the FOCUS example dataset D using defaults +#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE) +#' # Use starting parameters from parent only FOMC fit +#' fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, quiet = TRUE) +#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE, +#' parms.ini = fit.FOMC$bparms.ode) +#' +#' # Use stepwise fitting, using optimised parameters from parent only fit, SFORB +#' SFORB_SFO <- mkinmod( +#' parent = list(type = "SFORB", to = "m1", sink = TRUE), +#' m1 = list(type = "SFO")) +#' # Fit the model to the FOCUS example dataset D using defaults +#' fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, quiet = TRUE) +#' fit.SFORB_SFO.deSolve <- mkinfit(SFORB_SFO, FOCUS_2006_D, solution_type = "deSolve", +#' quiet = TRUE) +#' # Use starting parameters from parent only SFORB fit (not really needed in this case) +#' fit.SFORB = mkinfit("SFORB", FOCUS_2006_D, quiet = TRUE) +#' fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode, quiet = TRUE) +#' } +#' +#' \dontrun{ +#' # Weighted fits, including IRLS +#' SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), +#' m1 = mkinsub("SFO"), use_of_ff = "max") +#' f.noweight <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE) +#' summary(f.noweight) +#' f.obs <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "obs", quiet = TRUE) +#' summary(f.obs) +#' f.tc <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "tc", quiet = TRUE) +#' summary(f.tc) +#' } +#' +#' +#' @export +mkinfit <- function(mkinmod, observed, + parms.ini = "auto", + state.ini = "auto", + err.ini = "auto", + fixed_parms = NULL, + fixed_initials = names(mkinmod$diffs)[-1], + from_max_mean = FALSE, + solution_type = c("auto", "analytical", "eigen", "deSolve"), + method.ode = "lsoda", + use_compiled = "auto", + control = list(eval.max = 300, iter.max = 200), + transform_rates = TRUE, + transform_fractions = TRUE, + quiet = FALSE, + atol = 1e-8, rtol = 1e-10, n.outtimes = 100, + error_model = c("const", "obs", "tc"), + error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", "OLS"), + reweight.tol = 1e-8, reweight.max.iter = 10, + trace_parms = FALSE, + ...) +{ + # Check mkinmod and generate a model for the variable whith the highest value + # if a suitable string is given + parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic") + if (class(mkinmod) != "mkinmod") { + presumed_parent_name = observed[which.max(observed$value), "name"] + if (mkinmod[[1]] %in% parent_models_available) { + speclist <- list(list(type = mkinmod, sink = TRUE)) + names(speclist) <- presumed_parent_name + mkinmod <- mkinmod(speclist = speclist) + } else { + stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ", + paste(parent_models_available, collapse = ", ")) + } + } + + # Get the names of the state variables in the model + mod_vars <- names(mkinmod$diffs) + + # Get the names of observed variables + obs_vars <- names(mkinmod$spec) + + # Subset observed data with names of observed data in the model and remove NA values + observed <- subset(observed, name %in% obs_vars) + observed <- subset(observed, !is.na(value)) + + # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model) + if (any(observed$value == 0)) { + warning("Observations with value of zero were removed from the data") + observed <- subset(observed, value != 0) + } + + # Obtain data for decline from maximum mean value if requested + if (from_max_mean) { + # This is only used for simple decline models + if (length(obs_vars) > 1) + stop("Decline from maximum is only implemented for models with a single observed variable") + + means <- aggregate(value ~ time, data = observed, mean, na.rm=TRUE) + t_of_max <- means[which.max(means$value), "time"] + observed <- subset(observed, time >= t_of_max) + observed$time <- observed$time - t_of_max + } + + # Number observations used for fitting + n_observed <- nrow(observed) + + # Define starting values for parameters where not specified by the user + if (parms.ini[[1]] == "auto") parms.ini = vector() + + # Warn for inital parameter specifications that are not in the model + wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms) + if (length(wrongpar.names) > 0) { + warning("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "), + " not used in the model") + parms.ini <- parms.ini[setdiff(names(parms.ini), wrongpar.names)] + } + + # Warn that the sum of formation fractions may exceed one if they are not + # fitted in the transformed way + if (mkinmod$use_of_ff == "max" & transform_fractions == FALSE) { + warning("The sum of formation fractions may exceed one if you do not use ", + "transform_fractions = TRUE." ) + for (box in mod_vars) { + # Stop if formation fractions are not transformed and we have no sink + if (mkinmod$spec[[box]]$sink == FALSE) { + stop("If formation fractions are not transformed during the fitting, ", + "it is not supported to turn off pathways to sink.\n ", + "Consider turning on the transformation of formation fractions or ", + "setting up a model with use_of_ff = 'min'.\n") + } + } + } + + # Do not allow fixing formation fractions if we are using the ilr transformation, + # this is not supported + if (transform_fractions == TRUE && length(fixed_parms) > 0) { + if (any(grepl("^f_", fixed_parms))) { + stop("Fixing formation fractions is not supported when using the ilr ", + "transformation.") + } + } + + # Set initial parameter values, including a small increment (salt) + # to avoid linear dependencies (singular matrix) in Eigenvalue based solutions + k_salt = 0 + defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini)) + for (parmname in defaultpar.names) { + # Default values for rate constants, depending on the parameterisation + if (grepl("^k", parmname)) { + parms.ini[parmname] = 0.1 + k_salt + k_salt = k_salt + 1e-4 + } + # Default values for rate constants for reversible binding + if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1 + if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02 + # Default values for IORE exponents + if (grepl("^N", parmname)) parms.ini[parmname] = 1.1 + # Default values for the FOMC, DFOP and HS models + if (parmname == "alpha") parms.ini[parmname] = 1 + if (parmname == "beta") parms.ini[parmname] = 10 + if (parmname == "k1") parms.ini[parmname] = 0.1 + if (parmname == "k2") parms.ini[parmname] = 0.01 + if (parmname == "tb") parms.ini[parmname] = 5 + if (parmname == "g") parms.ini[parmname] = 0.5 + if (parmname == "kmax") parms.ini[parmname] = 0.1 + if (parmname == "k0") parms.ini[parmname] = 0.0001 + if (parmname == "r") parms.ini[parmname] = 0.2 + } + # Default values for formation fractions in case they are present + for (box in mod_vars) { + f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)] + if (length(f_names) > 0) { + # We need to differentiate between default and specified fractions + # and set the unspecified to 1 - sum(specified)/n_unspecified + f_default_names <- intersect(f_names, defaultpar.names) + f_specified_names <- setdiff(f_names, defaultpar.names) + sum_f_specified = sum(parms.ini[f_specified_names]) + if (sum_f_specified > 1) { + stop("Starting values for the formation fractions originating from ", + box, " sum up to more than 1.") + } + if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1 + else { + n_unspecified = length(f_default_names) + } + parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified + } + } + + # Set default for state.ini if appropriate + parent_name = names(mkinmod$spec)[[1]] + if (state.ini[1] == "auto") { + parent_time_0 = subset(observed, time == 0 & name == parent_name)$value + parent_time_0_mean = mean(parent_time_0, na.rm = TRUE) + if (is.na(parent_time_0_mean)) { + state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)) + } else { + state.ini = c(parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1)) + } + } + + # Name the inital state variable values if they are not named yet + if(is.null(names(state.ini))) names(state.ini) <- mod_vars + + # Transform initial parameter values for fitting + transparms.ini <- transform_odeparms(parms.ini, mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) + + # Parameters to be optimised: + # Kinetic parameters in parms.ini whose names are not in fixed_parms + parms.fixed <- parms.ini[fixed_parms] + parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)] + + transparms.fixed <- transform_odeparms(parms.fixed, mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) + transparms.optim <- transform_odeparms(parms.optim, mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) + + # Inital state variables in state.ini whose names are not in fixed_initials + state.ini.fixed <- state.ini[fixed_initials] + state.ini.optim <- state.ini[setdiff(names(state.ini), fixed_initials)] + + # Preserve names of state variables before renaming initial state variable + # parameters + state.ini.optim.boxnames <- names(state.ini.optim) + state.ini.fixed.boxnames <- names(state.ini.fixed) + if(length(state.ini.optim) > 0) { + names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_") + } + if(length(state.ini.fixed) > 0) { + names(state.ini.fixed) <- paste(names(state.ini.fixed), "0", sep="_") + } + + # Decide if the solution of the model can be based on a simple analytical + # formula, the spectral decomposition of the matrix (fundamental system) + # or a numeric ode solver from the deSolve package + # Prefer deSolve over eigen if a compiled model is present and use_compiled + # is not set to FALSE + solution_type = match.arg(solution_type) + if (solution_type == "analytical" && length(mkinmod$spec) > 1) + stop("Analytical solution not implemented for models with metabolites.") + if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat)) + stop("Eigenvalue based solution not possible, coefficient matrix not present.") + if (solution_type == "auto") { + if (length(mkinmod$spec) == 1) { + solution_type = "analytical" + } else { + if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) { + solution_type = "deSolve" + } else { + if (is.matrix(mkinmod$coefmat)) { + solution_type = "eigen" + if (max(observed$value, na.rm = TRUE) < 0.1) { + stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone") + } + } else { + solution_type = "deSolve" + } + } + } + } + + # Get the error model and the algorithm for fitting + err_mod <- match.arg(error_model) + error_model_algorithm = match.arg(error_model_algorithm) + if (error_model_algorithm == "OLS") { + if (err_mod != "const") stop("OLS is only appropriate for constant variance") + } + if (error_model_algorithm == "auto") { + error_model_algorithm = switch(err_mod, + const = "OLS", obs = "d_3", tc = "d_3") + } + errparm_names <- switch(err_mod, + "const" = "sigma", + "obs" = paste0("sigma_", obs_vars), + "tc" = c("sigma_low", "rsd_high")) + errparm_names_optim <- if (error_model_algorithm == "OLS") NULL else errparm_names + + # Define starting values for the error model + if (err.ini[1] != "auto") { + if (!identical(names(err.ini), errparm_names)) { + stop("Please supply initial values for error model components ", paste(errparm_names, collapse = ", ")) + } else { + errparms = err.ini + } + } else { + if (err_mod == "const") { + errparms = 3 + } + if (err_mod == "obs") { + errparms = rep(3, length(obs_vars)) + } + if (err_mod == "tc") { + errparms <- c(sigma_low = 0.1, rsd_high = 0.1) + } + names(errparms) <- errparm_names + } + if (error_model_algorithm == "OLS") { + errparms_optim <- NULL + } else { + errparms_optim <- errparms + } + + # Define outtimes for model solution. + # Include time points at which observed data are available + outtimes = sort(unique(c(observed$time, seq(min(observed$time), + max(observed$time), + length.out = n.outtimes)))) + + # Define the objective function for optimisation, including (back)transformations + cost_function <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...) + { + assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter + + # Trace parameter values if requested and if we are actually optimising + if(trace_parms & update_data) cat(P, "\n") + + # Determine local parameter values for the cost estimation + if (is.numeric(fixed_degparms)) { + cost_degparms <- fixed_degparms + cost_errparms <- P + degparms_fixed = TRUE + } else { + degparms_fixed = FALSE + } + + if (is.numeric(fixed_errparms)) { + cost_degparms <- P + cost_errparms <- fixed_errparms + errparms_fixed = TRUE + } else { + errparms_fixed = FALSE + } + + if (OLS) { + cost_degparms <- P + cost_errparms <- numeric(0) + } + + if (!OLS & !degparms_fixed & !errparms_fixed) { + cost_degparms <- P[1:(length(P) - length(errparms))] + cost_errparms <- P[(length(cost_degparms) + 1):length(P)] + } + + # Initial states for t0 + if(length(state.ini.optim) > 0) { + odeini <- c(cost_degparms[1:length(state.ini.optim)], state.ini.fixed) + names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames) + } else { + odeini <- state.ini.fixed + names(odeini) <- state.ini.fixed.boxnames + } + + odeparms.optim <- cost_degparms[(length(state.ini.optim) + 1):length(cost_degparms)] + + if (trans == TRUE) { + odeparms <- c(odeparms.optim, transparms.fixed) + parms <- backtransform_odeparms(odeparms, mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) + } else { + parms <- c(odeparms.optim, parms.fixed) + } + + # Solve the system with current parameter values + out <- mkinpredict(mkinmod, parms, + odeini, outtimes, + solution_type = solution_type, + use_compiled = use_compiled, + method.ode = method.ode, + atol = atol, rtol = rtol, ...) + + out_long <- mkin_wide_to_long(out, time = "time") + + if (err_mod == "const") { + observed$std <- if (OLS) NA else cost_errparms["sigma"] + } + if (err_mod == "obs") { + std_names <- paste0("sigma_", observed$name) + observed$std <- cost_errparms[std_names] + } + if (err_mod == "tc") { + tmp <- merge(observed, out_long, by = c("time", "name")) + tmp$name <- ordered(tmp$name, levels = obs_vars) + tmp <- tmp[order(tmp$name, tmp$time), ] + observed$std <- sqrt(cost_errparms["sigma_low"]^2 + tmp$value.y^2 * cost_errparms["rsd_high"]^2) + } + + cost_data <- merge(observed[c("name", "time", "value", "std")], out_long, + by = c("name", "time"), suffixes = c(".observed", ".predicted")) + + if (OLS) { + # Cost is the sum of squared residuals + cost <- with(cost_data, sum((value.observed - value.predicted)^2)) + } else { + # Cost is the negative log-likelihood + cost <- - with(cost_data, + sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE))) + } + + # We update the current cost and data during the optimisation, not + # during hessian calculations + if (update_data) { + + assign("out_predicted", out_long, inherits = TRUE) + assign("current_data", cost_data, inherits = TRUE) + + if (cost < cost.current) { + assign("cost.current", cost, inherits = TRUE) + if (!quiet) cat(ifelse(OLS, "Sum of squared residuals", "Negative log-likelihood"), + " at call ", calls, ": ", cost.current, "\n", sep = "") + } + } + return(cost) + } + + names_optim <- c(names(state.ini.optim), + names(transparms.optim), + errparm_names_optim) + n_optim <- length(names_optim) + + # Define lower and upper bounds other than -Inf and Inf for parameters + # for which no internal transformation is requested in the call to mkinfit + # and for optimised error model parameters + lower <- rep(-Inf, n_optim) + upper <- rep(Inf, n_optim) + names(lower) <- names(upper) <- names_optim + + # IORE exponents are not transformed, but need a lower bound + index_N <- grep("^N", names(lower)) + lower[index_N] <- 0 + + if (!transform_rates) { + index_k <- grep("^k_", names(lower)) + lower[index_k] <- 0 + index_k__iore <- grep("^k__iore_", names(lower)) + lower[index_k__iore] <- 0 + other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb", "r"), names(lower)) + lower[other_rate_parms] <- 0 + } + + if (!transform_fractions) { + index_f <- grep("^f_", names(upper)) + lower[index_f] <- 0 + upper[index_f] <- 1 + other_fraction_parms <- intersect(c("g"), names(upper)) + lower[other_fraction_parms] <- 0 + upper[other_fraction_parms] <- 1 + } + + if (err_mod == "const") { + if (error_model_algorithm != "OLS") { + lower["sigma"] <- 0 + } + } + if (err_mod == "obs") { + index_sigma <- grep("^sigma_", names(lower)) + lower[index_sigma] <- 0 + } + if (err_mod == "tc") { + lower["sigma_low"] <- 0 + lower["rsd_high"] <- 0 + } + + # Counter for cost function evaluations + calls = 0 + cost.current <- Inf + out_predicted <- NA + current_data <- NA + + # Show parameter names if tracing is requested + if(trace_parms) cat(names_optim, "\n") + + # browser() + + # Do the fit and take the time until the hessians are calculated + fit_time <- system.time({ + degparms <- c(state.ini.optim, transparms.optim) + n_degparms <- length(degparms) + degparms_index <- seq(1, n_degparms) + errparms_index <- seq(n_degparms + 1, length.out = length(errparms)) + + if (error_model_algorithm == "d_3") { + if (!quiet) message("Directly optimising the complete model") + parms.start <- c(degparms, errparms) + fit_direct <- nlminb(parms.start, cost_function, + lower = lower[names(parms.start)], + upper = upper[names(parms.start)], + control = control, ...) + fit_direct$logLik <- - cost.current + if (error_model_algorithm == "direct") { + degparms <- fit_direct$par[degparms_index] + errparms <- fit_direct$par[errparms_index] + } else { + cost.current <- Inf # reset to avoid conflict with the OLS step + } + } + if (error_model_algorithm != "direct") { + if (!quiet) message("Ordinary least squares optimisation") + fit <- nlminb(degparms, cost_function, control = control, + lower = lower[names(degparms)], + upper = upper[names(degparms)], OLS = TRUE, ...) + degparms <- fit$par + + # Get the maximum likelihood estimate for sigma at the optimum parameter values + current_data$residual <- current_data$value.observed - current_data$value.predicted + sigma_mle <- sqrt(sum(current_data$residual^2)/nrow(current_data)) + + # Use that estimate for the constant variance, or as first guess if err_mod = "obs" + if (err_mod != "tc") { + errparms[names(errparms)] <- sigma_mle + } + fit$par <- c(fit$par, errparms) + + cost.current <- cost_function(c(degparms, errparms), OLS = FALSE) + fit$logLik <- - cost.current + } + if (error_model_algorithm %in% c("threestep", "fourstep", "d_3")) { + if (!quiet) message("Optimising the error model") + fit <- nlminb(errparms, cost_function, control = control, + lower = lower[names(errparms)], + upper = upper[names(errparms)], + fixed_degparms = degparms, ...) + errparms <- fit$par + } + if (error_model_algorithm == "fourstep") { + if (!quiet) message("Optimising the degradation model") + fit <- nlminb(degparms, cost_function, control = control, + lower = lower[names(degparms)], + upper = upper[names(degparms)], + fixed_errparms = errparms, ...) + degparms <- fit$par + } + if (error_model_algorithm %in% + c("direct", "twostep", "threestep", "fourstep", "d_3")) { + if (!quiet) message("Optimising the complete model") + parms.start <- c(degparms, errparms) + fit <- nlminb(parms.start, cost_function, + lower = lower[names(parms.start)], + upper = upper[names(parms.start)], + control = control, ...) + degparms <- fit$par[degparms_index] + errparms <- fit$par[errparms_index] + fit$logLik <- - cost.current + + if (error_model_algorithm == "d_3") { + d_3_messages = c( + same = "Direct fitting and three-step fitting yield approximately the same likelihood", + threestep = "Three-step fitting yielded a higher likelihood than direct fitting", + direct = "Direct fitting yielded a higher likelihood than three-step fitting") + rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik)) + if (rel_diff < 0.0001) { + if (!quiet) message(d_3_messages["same"]) + fit$d_3_message <- d_3_messages["same"] + } else { + if (fit$logLik > fit_direct$logLik) { + if (!quiet) message(d_3_messages["threestep"]) + fit$d_3_message <- d_3_messages["threestep"] + } else { + if (!quiet) message(d_3_messages["direct"]) + fit <- fit_direct + fit$d_3_message <- d_3_messages["direct"] + } + } + } + } + if (err_mod != "const" & error_model_algorithm == "IRLS") { + reweight.diff <- 1 + n.iter <- 0 + errparms_last <- errparms + + while (reweight.diff > reweight.tol & + n.iter < reweight.max.iter) { + + if (!quiet) message("Optimising the error model") + fit <- nlminb(errparms, cost_function, control = control, + lower = lower[names(errparms)], + upper = upper[names(errparms)], + fixed_degparms = degparms, ...) + errparms <- fit$par + + if (!quiet) message("Optimising the degradation model") + fit <- nlminb(degparms, cost_function, control = control, + lower = lower[names(degparms)], + upper = upper[names(degparms)], + fixed_errparms = errparms, ...) + degparms <- fit$par + + reweight.diff <- dist(rbind(errparms, errparms_last)) + errparms_last <- errparms + + fit$par <- c(fit$par, errparms) + cost.current <- cost_function(c(degparms, errparms), OLS = FALSE) + fit$logLik <- - cost.current + } + } + + fit$hessian <- try(numDeriv::hessian(cost_function, c(degparms, errparms), OLS = FALSE, + update_data = FALSE), silent = TRUE) + + # Backtransform parameters + bparms.optim = backtransform_odeparms(fit$par, mkinmod, + transform_rates = transform_rates, + transform_fractions = transform_fractions) + bparms.fixed = c(state.ini.fixed, parms.fixed) + bparms.all = c(bparms.optim, parms.fixed) + + fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.all, errparms), + OLS = FALSE, trans = FALSE, update_data = FALSE), silent = TRUE) + }) + + fit$error_model_algorithm <- error_model_algorithm + + if (fit$convergence != 0) { + fit$warning = paste0("Optimisation did not converge:\n", fit$message) + warning(fit$warning) + } else { + if(!quiet) message("Optimisation successfully terminated.\n") + } + + # We need to return some more data for summary and plotting + fit$solution_type <- solution_type + fit$transform_rates <- transform_rates + fit$transform_fractions <- transform_fractions + fit$reweight.tol <- reweight.tol + fit$reweight.max.iter <- reweight.max.iter + fit$control <- control + fit$calls <- calls + fit$time <- fit_time + + # We also need the model for summary and plotting + fit$mkinmod <- mkinmod + + # We need data and predictions for summary and plotting + fit$observed <- observed + fit$obs_vars <- obs_vars + fit$predicted <- out_predicted + + # Residual sum of squares as a function of the fitted parameters + fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE) + + # Log-likelihood with possibility to fix degparms or errparms + fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE) { + - cost_function(P, fixed_degparms = fixed_degparms, + fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE) + } + + # Collect initial parameter values in three dataframes + fit$start <- data.frame(value = c(state.ini.optim, + parms.optim, errparms_optim)) + fit$start$type = c(rep("state", length(state.ini.optim)), + rep("deparm", length(parms.optim)), + rep("error", length(errparms_optim))) + + fit$start_transformed = data.frame( + value = c(state.ini.optim, transparms.optim, errparms_optim), + lower = lower, + upper = upper) + + fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed)) + fit$fixed$type = c(rep("state", length(state.ini.fixed)), + rep("deparm", length(parms.fixed))) + + # Sort observed, predicted and residuals + current_data$name <- ordered(current_data$name, levels = obs_vars) + + ordered_data <- current_data[order(current_data$name, current_data$time), ] + + fit$data <- data.frame(time = ordered_data$time, + variable = ordered_data$name, + observed = ordered_data$value.observed, + predicted = ordered_data$value.predicted) + + fit$data$residual <- fit$data$observed - fit$data$predicted + + fit$atol <- atol + fit$rtol <- rtol + fit$err_mod <- err_mod + + # Return different sets of backtransformed parameters for summary and plotting + fit$bparms.optim <- bparms.optim + fit$bparms.fixed <- bparms.fixed + + # Return ode and state parameters for further fitting + fit$bparms.ode <- bparms.all[mkinmod$parms] + fit$bparms.state <- c(bparms.all[setdiff(names(bparms.all), names(fit$bparms.ode))], + state.ini.fixed) + names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state)) + + fit$errparms <- errparms + fit$df.residual <- n_observed - length(c(degparms, errparms)) + + fit$date <- date() + fit$version <- as.character(utils::packageVersion("mkin")) + fit$Rversion <- paste(R.version$major, R.version$minor, sep=".") + + class(fit) <- c("mkinfit", "modFit") + return(fit) +} diff --git a/R/mkinmod.R b/R/mkinmod.R index 26148f18..a1ae0021 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,383 +1,483 @@ -# Copyright (C) 2010-2015,2019 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/> }}}
-
-mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verbose = FALSE)
-{
- if (is.null(speclist)) spec <- list(...)
- else spec <- speclist
- obs_vars <- names(spec)
-
- # Check if any of the names of the observed variables contains any other
- for (obs_var in obs_vars) {
- if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other")
- if (grepl("_to_", obs_var)) stop("Sorry, names of observed variables can not contain _to_")
- if (obs_var == "sink") stop("Naming a compound 'sink' is not supported")
- }
-
- if (!use_of_ff %in% c("min", "max"))
- stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'")
-
- # The returned model will be a list of character vectors, containing {{{
- # differential equations (if supported), parameter names and a mapping from
- # model variables to observed variables. If possible, a matrix representation
- # of the differential equations is included
- # Compiling the functions from the C code generated below only works if the
- # implicit assumption about differential equations specified below
- # is satisfied
- parms <- vector()
- # }}}
-
- # Do not return a coefficient matrix mat when FOMC, IORE, DFOP, HS or logistic is used for the parent {{{
- if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS", "logistic")) {
- mat = FALSE
- } else mat = TRUE
- #}}}
-
- # Establish a list of differential equations as well as a map from observed {{{
- # compartments to differential equations
- diffs <- vector()
- map <- list()
- for (varname in obs_vars)
- {
- # Check the type component of the compartment specification {{{
- if(is.null(spec[[varname]]$type)) stop(
- "Every part of the model specification must be a list containing a type component")
- if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB", "logistic")) stop(
- "Available types are SFO, FOMC, IORE, DFOP, HS, SFORB and logistic only")
- if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "logistic") & match(varname, obs_vars) != 1) {
- stop(paste("Types FOMC, DFOP, HS and logistic are only implemented for the first compartment,",
- "which is assumed to be the source compartment"))
- }
- #}}}
- # New (sub)compartments (boxes) needed for the model type {{{
- new_boxes <- switch(spec[[varname]]$type,
- SFO = varname,
- FOMC = varname,
- IORE = varname,
- DFOP = varname,
- HS = varname,
- logistic = varname,
- SFORB = paste(varname, c("free", "bound"), sep = "_")
- )
- map[[varname]] <- new_boxes
- names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes)) #}}}
- # Start a new differential equation for each new box {{{
- new_diffs <- paste("d_", new_boxes, " =", sep = "")
- names(new_diffs) <- new_boxes
- diffs <- c(diffs, new_diffs) #}}}
- } #}}}
-
- # Create content of differential equations and build parameter list {{{
- for (varname in obs_vars)
- {
- # Get the name of the box(es) we are working on for the decline term(s)
- box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB
- # Turn on sink if this is not explicitly excluded by the user by
- # specifying sink=FALSE
- if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
- if(spec[[varname]]$type %in% c("SFO", "IORE", "SFORB")) { # {{{ Add decline term
- if (use_of_ff == "min") { # Minimum use of formation fractions
- if(spec[[varname]]$type == "IORE" && length(spec[[varname]]$to) > 0) {
- stop("Transformation reactions from compounds modelled with IORE\n",
- "are only supported with formation fractions (use_of_ff = 'max')")
- }
- if(spec[[varname]]$sink) {
- # If sink is required, add first-order/IORE sink term
- k_compound_sink <- paste("k", box_1, "sink", sep = "_")
- if(spec[[varname]]$type == "IORE") {
- k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_")
- }
- parms <- c(parms, k_compound_sink)
- decline_term <- paste(k_compound_sink, "*", box_1)
- if(spec[[varname]]$type == "IORE") {
- N <- paste("N", box_1, sep = "_")
- parms <- c(parms, N)
- decline_term <- paste0(decline_term, "^", N)
- }
- } else { # otherwise no decline term needed here
- decline_term = "0"
- }
- } else {
- k_compound <- paste("k", box_1, sep = "_")
- if(spec[[varname]]$type == "IORE") {
- k_compound <- paste("k__iore", box_1, sep = "_")
- }
- parms <- c(parms, k_compound)
- decline_term <- paste(k_compound, "*", box_1)
- if(spec[[varname]]$type == "IORE") {
- N <- paste("N", box_1, sep = "_")
- parms <- c(parms, N)
- decline_term <- paste0(decline_term, "^", N)
- }
- }
- } #}}}
- if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term
- # From p. 53 of the FOCUS kinetics report, without the power function so it works in C
- decline_term <- paste("(alpha/beta) * 1/((time/beta) + 1) *", box_1)
- parms <- c(parms, "alpha", "beta")
- } #}}}
- if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term
- # From p. 57 of the FOCUS kinetics report
- decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1)
- parms <- c(parms, "k1", "k2", "g")
- } #}}}
- HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C
- if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term
- # From p. 55 of the FOCUS kinetics report
- decline_term <- paste(HS_decline, "*", box_1)
- parms <- c(parms, "k1", "k2", "tb")
- } #}}}
- if(spec[[varname]]$type == "logistic") { # {{{ Add logistic decline term
- # From p. 67 of the FOCUS kinetics report (2014)
- decline_term <- paste("(k0 * kmax)/(k0 + (kmax - k0) * exp(-r * time)) *", box_1)
- parms <- c(parms, "kmax", "k0", "r")
- } #}}}
- # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{
- diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}}
- if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms
- box_2 = map[[varname]][[2]]
- if (use_of_ff == "min") { # Minimum use of formation fractions
- k_free_bound <- paste("k", varname, "free", "bound", sep = "_")
- k_bound_free <- paste("k", varname, "bound", "free", sep = "_")
- parms <- c(parms, k_free_bound, k_bound_free)
- reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+",
- k_bound_free, "*", box_2)
- reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-",
- k_bound_free, "*", box_2)
- } else { # Use formation fractions also for the free compartment
- stop("The maximum use of formation fractions is not supported for SFORB models")
- # The problems were: Calculation of dissipation times did not work in this case
- # and the coefficient matrix is not generated correctly by the code present
- # in this file in this case
- #f_free_bound <- paste("f", varname, "free", "bound", sep = "_")
- #k_bound_free <- paste("k", varname, "bound", "free", sep = "_")
- #parms <- c(parms, f_free_bound, k_bound_free)
- #reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2)
- #reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-",
- # k_bound_free, "*", box_2)
- }
- diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1)
- diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2)
- } #}}}
-
- # Transfer between compartments#{{{
- to <- spec[[varname]]$to
- if(!is.null(to)) {
- # Name of box from which transfer takes place
- origin_box <- box_1
-
- # Number of targets
- n_targets = length(to)
-
- # Add transfer terms to listed compartments
- for (target in to) {
- if (!target %in% obs_vars) stop("You did not specify a submodel for target variable ", target)
- target_box <- switch(spec[[target]]$type,
- SFO = target,
- IORE = target,
- SFORB = paste(target, "free", sep = "_"))
- if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO", "SFORB"))
- {
- k_from_to <- paste("k", origin_box, target_box, sep = "_")
- parms <- c(parms, k_from_to)
- diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
- k_from_to, "*", origin_box)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- k_from_to, "*", origin_box)
- } else {
- # Do not introduce a formation fraction if this is the only target
- if (spec[[origin_box]]$sink == FALSE && n_targets == 1) {
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- decline_term)
- } else {
- fraction_to_target = paste("f", origin_box, "to", target, sep = "_")
- parms <- c(parms, fraction_to_target)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
- fraction_to_target, "*", decline_term)
- }
- }
- }
- } #}}}
- } #}}}
-
- model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff)
-
- # Create coefficient matrix if appropriate#{{{
- if (mat) {
- boxes <- names(diffs)
- n <- length(boxes)
- m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes))
-
- if (use_of_ff == "min") { # {{{ Minimum use of formation fractions
- for (from in boxes) {
- for (to in boxes) {
- if (from == to) { # diagonal elements
- k.candidate = paste("k", from, c(boxes, "sink"), sep = "_")
- k.candidate = sub("free.*bound", "free_bound", k.candidate)
- k.candidate = sub("bound.*free", "bound_free", k.candidate)
- k.effective = intersect(model$parms, k.candidate)
- m[from,to] = ifelse(length(k.effective) > 0,
- paste("-", k.effective, collapse = " "), "0")
-
- } else { # off-diagonal elements
- k.candidate = paste("k", from, to, sep = "_")
- if (sub("_free$", "", from) == sub("_bound$", "", to)) {
- k.candidate = paste("k", sub("_free$", "_free_bound", from), sep = "_")
- }
- if (sub("_bound$", "", from) == sub("_free$", "", to)) {
- k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep = "_")
- }
- k.effective = intersect(model$parms, k.candidate)
- m[to, from] = ifelse(length(k.effective) > 0,
- k.effective, "0")
- }
- }
- } # }}}
- } else { # {{{ Use formation fractions where possible
- for (from in boxes) {
- for (to in boxes) {
- if (from == to) { # diagonal elements
- k.candidate = paste("k", from, sep = "_")
- m[from,to] = ifelse(k.candidate %in% model$parms,
- paste("-", k.candidate), "0")
- if(grepl("_free", from)) { # add transfer to bound compartment for SFORB
- m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep = "_"))
- }
- if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB
- m[from,to] = paste("- k", from, "free", sep = "_")
- }
- m[from,to] = m[from,to]
- } else { # off-diagonal elements
- f.candidate = paste("f", from, "to", to, sep = "_")
- k.candidate = paste("k", from, to, sep = "_")
- # SFORB with maximum use of formation fractions not implemented, see above
- m[to, from] = ifelse(f.candidate %in% model$parms,
- paste(f.candidate, " * k_", from, sep = ""),
- ifelse(k.candidate %in% model$parms, k.candidate, "0"))
- # Special case: singular pathway and no sink
- if (spec[[from]]$sink == FALSE && length(spec[[from]]$to) == 1 && to %in% spec[[from]]$to) {
- m[to, from] = paste("k", from, sep = "_")
- }
- }
- }
- }
- } # }}}
- model$coefmat <- m
- }#}}}
-
- # Try to create a function compiled from C code if more than one observed {{{
- # variable and gcc is available
- if (length(obs_vars) > 1) {
- if (Sys.which("gcc") != "") {
-
- # Translate the R code for the derivatives to C code
- diffs.C <- paste(diffs, collapse = ";\n")
- diffs.C <- paste0(diffs.C, ";")
-
- # HS
- diffs.C <- gsub(HS_decline, "(time <= tb ? k1 : k2)", diffs.C, fixed = TRUE)
-
- for (i in seq_along(diffs)) {
- state_var <- names(diffs)[i]
-
- # IORE
- if (state_var %in% obs_vars) {
- if (spec[[state_var]]$type == "IORE") {
- diffs.C <- gsub(paste0(state_var, "^N_", state_var),
- paste0("pow(y[", i - 1, "], N_", state_var, ")"),
- diffs.C, fixed = TRUE)
- }
- }
-
- # Replace d_... terms by f[i-1]
- # First line
- pattern <- paste0("^d_", state_var)
- replacement <- paste0("\nf[", i - 1, "]")
- diffs.C <- gsub(pattern, replacement, diffs.C)
- # Other lines
- pattern <- paste0("\\nd_", state_var)
- replacement <- paste0("\nf[", i - 1, "]")
- diffs.C <- gsub(pattern, replacement, diffs.C)
-
- # Replace names of observed variables by y[i],
- # making the implicit assumption that the observed variables only occur after "* "
- pattern <- paste0("\\* ", state_var)
- replacement <- paste0("* y[", i - 1, "]")
- diffs.C <- gsub(pattern, replacement, diffs.C)
- }
-
- derivs_sig <- signature(n = "integer", t = "numeric", y = "numeric",
- f = "numeric", rpar = "numeric", ipar = "integer")
-
- # Declare the time variable in the body of the function if it is used
- derivs_code <- if (spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) {
- paste0("double time = *t;\n", diffs.C)
- } else {
- diffs.C
- }
-
- # Define the function initializing the parameters
- npar <- length(parms)
- initpar_code <- paste0(
- "static double parms [", npar, "];\n",
- paste0("#define ", parms, " parms[", 0:(npar - 1), "]\n", collapse = ""),
- "\n",
- "void initpar(void (* odeparms)(int *, double *)) {\n",
- " int N = ", npar, ";\n",
- " odeparms(&N, parms);\n",
- "}\n\n")
-
- # Try to build a shared library
- cf <- try(cfunction(list(func = derivs_sig), derivs_code,
- otherdefs = initpar_code,
- verbose = verbose,
- convention = ".C", language = "C"),
- silent = TRUE)
-
- if (!inherits(cf, "try-error")) {
- if (!quiet) message("Successfully compiled differential equation model from auto-generated C code.")
- model$cf <- cf
- }
- }
- }
- # }}}
-
- class(model) <- "mkinmod"
- return(model)
-}
-
-print.mkinmod <- function(x, ...) {
- cat("<mkinmod> model generated with\n")
- cat("Use of formation fractions $use_of_ff:", x$use_of_ff, "\n")
- cat("Specification $spec:\n")
- for (obs in names(x$spec)) {
- cat("$", obs, "\n", sep = "")
- spl <- x$spec[[obs]]
- cat("$type:", spl$type)
- if (!is.null(spl$to) && length(spl$to)) cat("; $to: ", paste(spl$to, collapse = ", "), sep = "")
- cat("; $sink: ", spl$sink, sep = "")
- if (!is.null(spl$full_name)) if (!is.na(spl$full_name)) cat("; $full_name:", spl$full_name)
- cat("\n")
- }
- if (is.matrix(x$coefmat)) cat("Coefficient matrix $coefmat available\n")
- if (!is.null(x$cf)) cat("Compiled model $cf available\n")
- cat("Differential equations:\n")
- nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
- writeLines(strwrap(nice_diffs, exdent = 11))
-}
-# vim: set foldmethod=marker ts=2 sw=2 expandtab:
+#' Function to set up a kinetic model with one or more state variables +#' +#' The function usually takes several expressions, each assigning a compound +#' name to a list, specifying the kinetic model type and reaction or transfer +#' to other observed compartments. Instead of specifying several expressions, a +#' list of lists can be given in the speclist argument. +#' +#' For the definition of model types and their parameters, the equations given +#' in the FOCUS and NAFTA guidance documents are used. +#' +#' @param ... For each observed variable, a list has to be specified as an +#' argument, containing at least a component \code{type}, specifying the type +#' of kinetics to use for the variable. Currently, single first order +#' kinetics "SFO", indeterminate order rate equation kinetics "IORE", or +#' single first order with reversible binding "SFORB" are implemented for all +#' variables, while "FOMC", "DFOP" and "HS" can additionally be chosen for +#' the first variable which is assumed to be the source compartment. +#' Additionally, each component of the list can include a character vector +#' \code{to}, specifying names of variables to which a transfer is to be +#' assumed in the model. If the argument \code{use_of_ff} is set to "min" +#' (default) and the model for the compartment is "SFO" or "SFORB", an +#' additional component of the list can be "sink=FALSE" effectively fixing +#' the flux to sink to zero. +#' @param use_of_ff Specification of the use of formation fractions in the +#' model equations and, if applicable, the coefficient matrix. If "min", a +#' minimum use of formation fractions is made in order to avoid fitting the +#' product of formation fractions and rate constants. If "max", formation +#' fractions are always used. +#' @param speclist The specification of the observed variables and their +#' submodel types and pathways can be given as a single list using this +#' argument. Default is NULL. +#' @param quiet Should messages be suppressed? +#' @param verbose If \code{TRUE}, passed to \code{\link{cfunction}} if +#' applicable to give detailed information about the C function being built. +#' @importFrom methods signature +#' @importFrom inline cfunction +#' @return A list of class \code{mkinmod} for use with \code{\link{mkinfit}}, +#' containing, among others, +#' \item{diffs}{ +#' A vector of string representations of differential equations, one for +#' each modelling variable. +#' } +#' \item{map}{ +#' A list containing named character vectors for each observed variable, +#' specifying the modelling variables by which it is represented. +#' } +#' \item{use_of_ff}{ +#' The content of \code{use_of_ff} is passed on in this list component. +#' } +#' \item{coefmat}{ +#' The coefficient matrix, if the system of differential equations can be +#' represented by one. +#' } +#' \item{ll}{ +#' The likelihood function, taking the parameter vector as the first argument. +#' } +#' @note The IORE submodel is not well tested for metabolites. When using this +#' model for metabolites, you may want to read the second note in the help +#' page to \code{\link{mkinfit}}. +#' @author Johannes Ranke +#' @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} +#' +#' NAFTA Technical Working Group on Pesticides (not dated) Guidance for +#' Evaluating and Calculating Degradation Kinetics in Environmental Media +#' @examples +#' +#' # Specify the SFO model (this is not needed any more, as we can now mkinfit("SFO", ...) +#' SFO <- mkinmod(parent = list(type = "SFO")) +#' +#' # One parent compound, one metabolite, both single first order +#' SFO_SFO <- mkinmod( +#' parent = mkinsub("SFO", "m1"), +#' m1 = mkinsub("SFO")) +#' +#' \dontrun{ +#' # The above model used to be specified like this, before the advent of mkinsub() +#' SFO_SFO <- mkinmod( +#' parent = list(type = "SFO", to = "m1"), +#' m1 = list(type = "SFO")) +#' +#' # Show details of creating the C function +#' SFO_SFO <- mkinmod( +#' parent = mkinsub("SFO", "m1"), +#' m1 = mkinsub("SFO"), verbose = TRUE) +#' +#' # If we have several parallel metabolites +#' # (compare tests/testthat/test_synthetic_data_for_UBA_2014.R) +#' m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")), +#' M1 = mkinsub("SFO"), +#' M2 = mkinsub("SFO"), +#' use_of_ff = "max", quiet = TRUE) +#' +#' fit_DFOP_par_c <- mkinfit(m_synth_DFOP_par, +#' synthetic_data_for_UBA_2014[[12]]$data, +#' quiet = TRUE) +#' } +#' +#' @export mkinmod +mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verbose = FALSE) +{ + if (is.null(speclist)) spec <- list(...) + else spec <- speclist + obs_vars <- names(spec) + + # Check if any of the names of the observed variables contains any other + for (obs_var in obs_vars) { + if (length(grep(obs_var, obs_vars)) > 1) stop("Sorry, variable names can not contain each other") + if (grepl("_to_", obs_var)) stop("Sorry, names of observed variables can not contain _to_") + if (obs_var == "sink") stop("Naming a compound 'sink' is not supported") + } + + if (!use_of_ff %in% c("min", "max")) + stop("The use of formation fractions 'use_of_ff' can only be 'min' or 'max'") + + # The returned model will be a list of character vectors, containing {{{ + # differential equations (if supported), parameter names and a mapping from + # model variables to observed variables. If possible, a matrix representation + # of the differential equations is included + # Compiling the functions from the C code generated below only works if the + # implicit assumption about differential equations specified below + # is satisfied + parms <- vector() + # }}} + + # Do not return a coefficient matrix mat when FOMC, IORE, DFOP, HS or logistic is used for the parent {{{ + if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS", "logistic")) { + mat = FALSE + } else mat = TRUE + #}}} + + # Establish a list of differential equations as well as a map from observed {{{ + # compartments to differential equations + diffs <- vector() + map <- list() + for (varname in obs_vars) + { + # Check the type component of the compartment specification {{{ + if(is.null(spec[[varname]]$type)) stop( + "Every part of the model specification must be a list containing a type component") + if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB", "logistic")) stop( + "Available types are SFO, FOMC, IORE, DFOP, HS, SFORB and logistic only") + if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS", "logistic") & match(varname, obs_vars) != 1) { + stop(paste("Types FOMC, DFOP, HS and logistic are only implemented for the first compartment,", + "which is assumed to be the source compartment")) + } + #}}} + # New (sub)compartments (boxes) needed for the model type {{{ + new_boxes <- switch(spec[[varname]]$type, + SFO = varname, + FOMC = varname, + IORE = varname, + DFOP = varname, + HS = varname, + logistic = varname, + SFORB = paste(varname, c("free", "bound"), sep = "_") + ) + map[[varname]] <- new_boxes + names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes)) #}}} + # Start a new differential equation for each new box {{{ + new_diffs <- paste("d_", new_boxes, " =", sep = "") + names(new_diffs) <- new_boxes + diffs <- c(diffs, new_diffs) #}}} + } #}}} + + # Create content of differential equations and build parameter list {{{ + for (varname in obs_vars) + { + # Get the name of the box(es) we are working on for the decline term(s) + box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB + # Turn on sink if this is not explicitly excluded by the user by + # specifying sink=FALSE + if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE + if(spec[[varname]]$type %in% c("SFO", "IORE", "SFORB")) { # {{{ Add decline term + if (use_of_ff == "min") { # Minimum use of formation fractions + if(spec[[varname]]$type == "IORE" && length(spec[[varname]]$to) > 0) { + stop("Transformation reactions from compounds modelled with IORE\n", + "are only supported with formation fractions (use_of_ff = 'max')") + } + if(spec[[varname]]$sink) { + # If sink is required, add first-order/IORE sink term + k_compound_sink <- paste("k", box_1, "sink", sep = "_") + if(spec[[varname]]$type == "IORE") { + k_compound_sink <- paste("k__iore", box_1, "sink", sep = "_") + } + parms <- c(parms, k_compound_sink) + decline_term <- paste(k_compound_sink, "*", box_1) + if(spec[[varname]]$type == "IORE") { + N <- paste("N", box_1, sep = "_") + parms <- c(parms, N) + decline_term <- paste0(decline_term, "^", N) + } + } else { # otherwise no decline term needed here + decline_term = "0" + } + } else { + k_compound <- paste("k", box_1, sep = "_") + if(spec[[varname]]$type == "IORE") { + k_compound <- paste("k__iore", box_1, sep = "_") + } + parms <- c(parms, k_compound) + decline_term <- paste(k_compound, "*", box_1) + if(spec[[varname]]$type == "IORE") { + N <- paste("N", box_1, sep = "_") + parms <- c(parms, N) + decline_term <- paste0(decline_term, "^", N) + } + } + } #}}} + if(spec[[varname]]$type == "FOMC") { # {{{ Add FOMC decline term + # From p. 53 of the FOCUS kinetics report, without the power function so it works in C + decline_term <- paste("(alpha/beta) * 1/((time/beta) + 1) *", box_1) + parms <- c(parms, "alpha", "beta") + } #}}} + if(spec[[varname]]$type == "DFOP") { # {{{ Add DFOP decline term + # From p. 57 of the FOCUS kinetics report + decline_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1) + parms <- c(parms, "k1", "k2", "g") + } #}}} + HS_decline <- "ifelse(time <= tb, k1, k2)" # Used below for automatic translation to C + if(spec[[varname]]$type == "HS") { # {{{ Add HS decline term + # From p. 55 of the FOCUS kinetics report + decline_term <- paste(HS_decline, "*", box_1) + parms <- c(parms, "k1", "k2", "tb") + } #}}} + if(spec[[varname]]$type == "logistic") { # {{{ Add logistic decline term + # From p. 67 of the FOCUS kinetics report (2014) + decline_term <- paste("(k0 * kmax)/(k0 + (kmax - k0) * exp(-r * time)) *", box_1) + parms <- c(parms, "kmax", "k0", "r") + } #}}} + # Add origin decline term to box 1 (usually the only box, unless type is SFORB)#{{{ + diffs[[box_1]] <- paste(diffs[[box_1]], "-", decline_term)#}}} + if(spec[[varname]]$type == "SFORB") { # {{{ Add SFORB reversible binding terms + box_2 = map[[varname]][[2]] + if (use_of_ff == "min") { # Minimum use of formation fractions + k_free_bound <- paste("k", varname, "free", "bound", sep = "_") + k_bound_free <- paste("k", varname, "bound", "free", sep = "_") + parms <- c(parms, k_free_bound, k_bound_free) + reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+", + k_bound_free, "*", box_2) + reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-", + k_bound_free, "*", box_2) + } else { # Use formation fractions also for the free compartment + stop("The maximum use of formation fractions is not supported for SFORB models") + # The problems were: Calculation of dissipation times did not work in this case + # and the coefficient matrix is not generated correctly by the code present + # in this file in this case + #f_free_bound <- paste("f", varname, "free", "bound", sep = "_") + #k_bound_free <- paste("k", varname, "bound", "free", sep = "_") + #parms <- c(parms, f_free_bound, k_bound_free) + #reversible_binding_term_1 <- paste("+", k_bound_free, "*", box_2) + #reversible_binding_term_2 <- paste("+", f_free_bound, "*", k_compound, "*", box_1, "-", + # k_bound_free, "*", box_2) + } + diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1) + diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2) + } #}}} + + # Transfer between compartments#{{{ + to <- spec[[varname]]$to + if(!is.null(to)) { + # Name of box from which transfer takes place + origin_box <- box_1 + + # Number of targets + n_targets = length(to) + + # Add transfer terms to listed compartments + for (target in to) { + if (!target %in% obs_vars) stop("You did not specify a submodel for target variable ", target) + target_box <- switch(spec[[target]]$type, + SFO = target, + IORE = target, + SFORB = paste(target, "free", sep = "_")) + if (use_of_ff == "min" && spec[[varname]]$type %in% c("SFO", "SFORB")) + { + k_from_to <- paste("k", origin_box, target_box, sep = "_") + parms <- c(parms, k_from_to) + diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", + k_from_to, "*", origin_box) + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + k_from_to, "*", origin_box) + } else { + # Do not introduce a formation fraction if this is the only target + if (spec[[origin_box]]$sink == FALSE && n_targets == 1) { + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + decline_term) + } else { + fraction_to_target = paste("f", origin_box, "to", target, sep = "_") + parms <- c(parms, fraction_to_target) + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + fraction_to_target, "*", decline_term) + } + } + } + } #}}} + } #}}} + + model <- list(diffs = diffs, parms = parms, map = map, spec = spec, use_of_ff = use_of_ff) + + # Create coefficient matrix if appropriate#{{{ + if (mat) { + boxes <- names(diffs) + n <- length(boxes) + m <- matrix(nrow=n, ncol=n, dimnames=list(boxes, boxes)) + + if (use_of_ff == "min") { # {{{ Minimum use of formation fractions + for (from in boxes) { + for (to in boxes) { + if (from == to) { # diagonal elements + k.candidate = paste("k", from, c(boxes, "sink"), sep = "_") + k.candidate = sub("free.*bound", "free_bound", k.candidate) + k.candidate = sub("bound.*free", "bound_free", k.candidate) + k.effective = intersect(model$parms, k.candidate) + m[from,to] = ifelse(length(k.effective) > 0, + paste("-", k.effective, collapse = " "), "0") + + } else { # off-diagonal elements + k.candidate = paste("k", from, to, sep = "_") + if (sub("_free$", "", from) == sub("_bound$", "", to)) { + k.candidate = paste("k", sub("_free$", "_free_bound", from), sep = "_") + } + if (sub("_bound$", "", from) == sub("_free$", "", to)) { + k.candidate = paste("k", sub("_bound$", "_bound_free", from), sep = "_") + } + k.effective = intersect(model$parms, k.candidate) + m[to, from] = ifelse(length(k.effective) > 0, + k.effective, "0") + } + } + } # }}} + } else { # {{{ Use formation fractions where possible + for (from in boxes) { + for (to in boxes) { + if (from == to) { # diagonal elements + k.candidate = paste("k", from, sep = "_") + m[from,to] = ifelse(k.candidate %in% model$parms, + paste("-", k.candidate), "0") + if(grepl("_free", from)) { # add transfer to bound compartment for SFORB + m[from,to] = paste(m[from,to], "-", paste("k", from, "bound", sep = "_")) + } + if(grepl("_bound", from)) { # add backtransfer to free compartment for SFORB + m[from,to] = paste("- k", from, "free", sep = "_") + } + m[from,to] = m[from,to] + } else { # off-diagonal elements + f.candidate = paste("f", from, "to", to, sep = "_") + k.candidate = paste("k", from, to, sep = "_") + # SFORB with maximum use of formation fractions not implemented, see above + m[to, from] = ifelse(f.candidate %in% model$parms, + paste(f.candidate, " * k_", from, sep = ""), + ifelse(k.candidate %in% model$parms, k.candidate, "0")) + # Special case: singular pathway and no sink + if (spec[[from]]$sink == FALSE && length(spec[[from]]$to) == 1 && to %in% spec[[from]]$to) { + m[to, from] = paste("k", from, sep = "_") + } + } + } + } + } # }}} + model$coefmat <- m + }#}}} + + # Try to create a function compiled from C code if more than one observed {{{ + # variable and gcc is available + if (length(obs_vars) > 1) { + if (Sys.which("gcc") != "") { + + # Translate the R code for the derivatives to C code + diffs.C <- paste(diffs, collapse = ";\n") + diffs.C <- paste0(diffs.C, ";") + + # HS + diffs.C <- gsub(HS_decline, "(time <= tb ? k1 : k2)", diffs.C, fixed = TRUE) + + for (i in seq_along(diffs)) { + state_var <- names(diffs)[i] + + # IORE + if (state_var %in% obs_vars) { + if (spec[[state_var]]$type == "IORE") { + diffs.C <- gsub(paste0(state_var, "^N_", state_var), + paste0("pow(y[", i - 1, "], N_", state_var, ")"), + diffs.C, fixed = TRUE) + } + } + + # Replace d_... terms by f[i-1] + # First line + pattern <- paste0("^d_", state_var) + replacement <- paste0("\nf[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + # Other lines + pattern <- paste0("\\nd_", state_var) + replacement <- paste0("\nf[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + + # Replace names of observed variables by y[i], + # making the implicit assumption that the observed variables only occur after "* " + pattern <- paste0("\\* ", state_var) + replacement <- paste0("* y[", i - 1, "]") + diffs.C <- gsub(pattern, replacement, diffs.C) + } + + derivs_sig <- signature(n = "integer", t = "numeric", y = "numeric", + f = "numeric", rpar = "numeric", ipar = "integer") + + # Declare the time variable in the body of the function if it is used + derivs_code <- if (spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) { + paste0("double time = *t;\n", diffs.C) + } else { + diffs.C + } + + # Define the function initializing the parameters + npar <- length(parms) + initpar_code <- paste0( + "static double parms [", npar, "];\n", + paste0("#define ", parms, " parms[", 0:(npar - 1), "]\n", collapse = ""), + "\n", + "void initpar(void (* odeparms)(int *, double *)) {\n", + " int N = ", npar, ";\n", + " odeparms(&N, parms);\n", + "}\n\n") + + # Try to build a shared library + cf <- try(cfunction(list(func = derivs_sig), derivs_code, + otherdefs = initpar_code, + verbose = verbose, + convention = ".C", language = "C"), + silent = TRUE) + + if (!inherits(cf, "try-error")) { + if (!quiet) message("Successfully compiled differential equation model from auto-generated C code.") + model$cf <- cf + } + } + } + # }}} + + class(model) <- "mkinmod" + return(model) +} + +#' Print mkinmod objects +#' +#' Print mkinmod objects in a way that the user finds his way to get to its +#' components. +#' +#' @param x An \code{\link{mkinmod}} object. +#' @param \dots Not used. +#' @examples +#' +#' m_synth_SFO_lin <- mkinmod(parent = list(type = "SFO", to = "M1"), +#' M1 = list(type = "SFO", to = "M2"), +#' M2 = list(type = "SFO"), use_of_ff = "max") +#' +#' print(m_synth_SFO_lin) +#' +#' @export +print.mkinmod <- function(x, ...) { + cat("<mkinmod> model generated with\n") + cat("Use of formation fractions $use_of_ff:", x$use_of_ff, "\n") + cat("Specification $spec:\n") + for (obs in names(x$spec)) { + cat("$", obs, "\n", sep = "") + spl <- x$spec[[obs]] + cat("$type:", spl$type) + if (!is.null(spl$to) && length(spl$to)) cat("; $to: ", paste(spl$to, collapse = ", "), sep = "") + cat("; $sink: ", spl$sink, sep = "") + if (!is.null(spl$full_name)) if (!is.na(spl$full_name)) cat("; $full_name:", spl$full_name) + cat("\n") + } + if (is.matrix(x$coefmat)) cat("Coefficient matrix $coefmat available\n") + if (!is.null(x$cf)) cat("Compiled model $cf available\n") + cat("Differential equations:\n") + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]]) + writeLines(strwrap(nice_diffs, exdent = 11)) +} +# vim: set foldmethod=marker ts=2 sw=2 expandtab: diff --git a/R/mkinparplot.R b/R/mkinparplot.R index af28e3a8..f9abab5b 100644 --- a/R/mkinparplot.R +++ b/R/mkinparplot.R @@ -1,20 +1,23 @@ -# Copyright (C) 2014 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/> +#' Function to plot the confidence intervals obtained using mkinfit +#' +#' This function plots the confidence intervals for the parameters fitted using +#' \code{\link{mkinfit}}. +#' +#' @param object A fit represented in an \code{\link{mkinfit}} object. +#' @return Nothing is returned by this function, as it is called for its side +#' effect, namely to produce a plot. +#' @author Johannes Ranke +#' @examples +#' +#' \dontrun{ +#' model <- mkinmod( +#' T245 = mkinsub("SFO", to = c("phenol"), sink = FALSE), +#' phenol = mkinsub("SFO", to = c("anisole")), +#' anisole = mkinsub("SFO"), use_of_ff = "max") +#' fit <- mkinfit(model, subset(mccall81_245T, soil == "Commerce"), quiet = TRUE) +#' mkinparplot(fit) +#' } +#' @export mkinparplot <- function(object) { state.optim = rownames(subset(object$start, type == "state")) deparms.optim = rownames(subset(object$start, type == "deparm")) diff --git a/R/mkinplot.R b/R/mkinplot.R deleted file mode 100644 index b9becfdf..00000000 --- a/R/mkinplot.R +++ /dev/null @@ -1,4 +0,0 @@ -mkinplot <- function(fit, ...) -{ - plot(fit, ...) -} diff --git a/R/mkinpredict.R b/R/mkinpredict.R index c36d724a..8949e800 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -1,22 +1,99 @@ -# Copyright (C) 2010-2016,2018,2019 Johannes Ranke -# Some lines in this code are copyright (C) 2013 Eurofins Regulatory AG -# 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/> - +#' Produce predictions from a kinetic model using specific parameters +#' +#' This function produces a time series for all the observed variables in a +#' kinetic model as specified by \code{\link{mkinmod}}, using a specific set of +#' kinetic parameters and initial values for the state variables. +#' +#' @aliases mkinpredict mkinpredict.mkinmod mkinpredict.mkinfit +#' @param x A kinetic model as produced by \code{\link{mkinmod}}, or a kinetic +#' fit as fitted by \code{\link{mkinfit}}. In the latter case, the fitted +#' parameters are used for the prediction. +#' @param odeparms A numeric vector specifying the parameters used in the +#' kinetic model, which is generally defined as a set of ordinary +#' differential equations. +#' @param odeini A numeric vectory containing the initial values of the state +#' variables of the model. Note that the state variables can differ from the +#' observed variables, for example in the case of the SFORB model. +#' @param outtimes A numeric vector specifying the time points for which model +#' predictions should be generated. +#' @param solution_type The method that should be used for producing the +#' predictions. This should generally be "analytical" if there is only one +#' observed variable, and usually "deSolve" in the case of several observed +#' variables. The third possibility "eigen" is faster but not applicable to +#' some models e.g. using FOMC for the parent compound. +#' @param method.ode The solution method passed via \code{\link{mkinpredict}} +#' to \code{\link{ode}} in case the solution type is "deSolve". The default +#' "lsoda" is performant, but sometimes fails to converge. +#' @param use_compiled If set to \code{FALSE}, no compiled version of the +#' \code{\link{mkinmod}} model is used, even if is present. +#' @param atol Absolute error tolerance, passed to \code{\link{ode}}. Default +#' is 1e-8, lower than in \code{\link{lsoda}}. +#' @param rtol Absolute error tolerance, passed to \code{\link{ode}}. Default +#' is 1e-10, much lower than in \code{\link{lsoda}}. +#' @param map_output Boolean to specify if the output should list values for +#' the observed variables (default) or for all state variables (if set to +#' FALSE). +#' @param \dots Further arguments passed to the ode solver in case such a +#' solver is used. +#' @import deSolve +#' @importFrom inline getDynLib +#' @return A matrix in the same format as the output of \code{\link{ode}}. +#' @author Johannes Ranke +#' @examples +#' +#' SFO <- mkinmod(degradinol = mkinsub("SFO")) +#' # Compare solution types +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "analytical") +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "deSolve") +#' mkinpredict(SFO, c(k_degradinol_sink = 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, +#' solution_type = "eigen") +#' +#' +#' # Compare integration methods to analytical solution +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' solution_type = "analytical")[21,] +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' method = "lsoda")[21,] +#' mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, +#' method = "ode45")[21,] +#' mkinpredict(SFO, c(k_degradinol_sink = 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), +#' seq(0, 20, by = 0.1))[201,] +#' mkinpredict(SFO, c(k_degradinol_sink = 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")) +#' system.time( +#' print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), +#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), +#' solution_type = "eigen")[201,])) +#' system.time( +#' print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), +#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), +#' solution_type = "deSolve")[201,])) +#' system.time( +#' print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), +#' c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), +#' solution_type = "deSolve", use_compiled = FALSE)[201,])) +#' +#' \dontrun{ +#' # Predict from a fitted model +#' f <- mkinfit(SFO_SFO, FOCUS_2006_C) +#' head(mkinpredict(f)) +#' } +#' +#' @export mkinpredict <- function(x, odeparms, odeini, outtimes = seq(0, 120, by = 0.1), solution_type = "deSolve", @@ -27,6 +104,8 @@ mkinpredict <- function(x, odeparms, odeini, UseMethod("mkinpredict", x) } +#' @rdname mkinpredict +#' @export mkinpredict.mkinmod <- function(x, odeparms = c(k_parent_sink = 0.1), odeini = c(parent = 100), @@ -164,6 +243,8 @@ mkinpredict.mkinmod <- function(x, } } +#' @rdname mkinpredict +#' @export mkinpredict.mkinfit <- function(x, odeparms = x$bparms.ode, odeini = x$bparms.state, diff --git a/R/mkinresplot.R b/R/mkinresplot.R index 974d0549..5377dbf2 100644 --- a/R/mkinresplot.R +++ b/R/mkinresplot.R @@ -1,65 +1,84 @@ -# Copyright (C) 2008-2014,2019 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/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual"))
-
-mkinresplot <- function (object,
- obs_vars = names(object$mkinmod$map),
- xlim = c(0, 1.1 * max(object$data$time)),
- xlab = "Time", ylab = "Residual",
- maxabs = "auto", legend= TRUE, lpos = "topright",
- col_obs = "auto", pch_obs = "auto",
- frame = TRUE,
- ...)
-{
- obs_vars_all <- as.character(unique(object$data$variable))
-
- if (length(obs_vars) > 0){
- obs_vars <- intersect(obs_vars_all, obs_vars)
- } else obs_vars <- obs_vars_all
-
- residuals <- subset(object$data, variable %in% obs_vars, residual)
-
- if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE)
-
- # Set colors and symbols
- if (col_obs[1] == "auto") {
- col_obs <- 1:length(obs_vars)
- }
-
- if (pch_obs[1] == "auto") {
- pch_obs <- 1:length(obs_vars)
- }
- names(col_obs) <- names(pch_obs) <- obs_vars
-
- plot(0, type = "n", frame = frame,
- xlab = xlab, ylab = ylab,
- xlim = xlim,
- ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...)
-
- for(obs_var in obs_vars){
- residuals_plot <- subset(object$data, variable == obs_var, c("time", "residual"))
- points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var])
- }
-
- abline(h = 0, lty = 2)
-
- if (legend == TRUE) {
- legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
- col = col_obs[obs_vars], pch = pch_obs[obs_vars])
- }
-}
+if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual")) + +#' Function to plot residuals stored in an mkin object +#' +#' This function plots the residuals for the specified subset of the observed +#' variables from an mkinfit object. A combined plot of the fitted model and +#' the residuals can be obtained using \code{\link{plot.mkinfit}} using the +#' argument \code{show_residuals = TRUE}. +#' +#' @param object A fit represented in an \code{\link{mkinfit}} object. +#' @param obs_vars A character vector of names of the observed variables for +#' which residuals should be plotted. Defaults to all observed variables in +#' the model +#' @param xlim plot range in x direction. +#' @param xlab Label for the x axis. Defaults to "Time [days]". +#' @param ylab Label for the y axis. Defaults to "Residual [\% of applied +#' radioactivity]". +#' @param maxabs Maximum absolute value of the residuals. This is used for the +#' scaling of the y axis and defaults to "auto". +#' @param legend Should a legend be plotted? Defaults to "TRUE". +#' @param lpos Where should the legend be placed? Default is "topright". Will +#' be passed on to \code{\link{legend}}. +#' @param col_obs Colors for the observed variables. +#' @param pch_obs Symbols to be used for the observed variables. +#' @param frame Should a frame be drawn around the plots? +#' @param \dots further arguments passed to \code{\link{plot}}. +#' @return Nothing is returned by this function, as it is called for its side +#' effect, namely to produce a plot. +#' @author Johannes Ranke +#' @seealso \code{\link{mkinplot}}, for a way to plot the data and the fitted +#' lines of the mkinfit object. +#' @examples +#' +#' model <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO")) +#' fit <- mkinfit(model, FOCUS_2006_D, quiet = TRUE) +#' mkinresplot(fit, "m1") +#' +#' @export +mkinresplot <- function (object, + obs_vars = names(object$mkinmod$map), + xlim = c(0, 1.1 * max(object$data$time)), + xlab = "Time", ylab = "Residual", + maxabs = "auto", legend= TRUE, lpos = "topright", + col_obs = "auto", pch_obs = "auto", + frame = TRUE, + ...) +{ + obs_vars_all <- as.character(unique(object$data$variable)) + + if (length(obs_vars) > 0){ + obs_vars <- intersect(obs_vars_all, obs_vars) + } else obs_vars <- obs_vars_all + + residuals <- subset(object$data, variable %in% obs_vars, residual) + + if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) + + # Set colors and symbols + if (col_obs[1] == "auto") { + col_obs <- 1:length(obs_vars) + } + + if (pch_obs[1] == "auto") { + pch_obs <- 1:length(obs_vars) + } + names(col_obs) <- names(pch_obs) <- obs_vars + + plot(0, type = "n", frame = frame, + xlab = xlab, ylab = ylab, + xlim = xlim, + ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...) + + for(obs_var in obs_vars){ + residuals_plot <- subset(object$data, variable == obs_var, c("time", "residual")) + points(residuals_plot, pch = pch_obs[obs_var], col = col_obs[obs_var]) + } + + abline(h = 0, lty = 2) + + if (legend == TRUE) { + legend(lpos, inset = c(0.05, 0.05), legend = obs_vars, + col = col_obs[obs_vars], pch = pch_obs[obs_vars]) + } +} diff --git a/R/mkinsub.R b/R/mkinsub.R index 99c3ea20..db91ca00 100644 --- a/R/mkinsub.R +++ b/R/mkinsub.R @@ -1,23 +1,39 @@ -# Copyright (C) 2014,2015 Johannes Ranke -# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG -# 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/> +#' Function to set up a kinetic submodel for one state variable +#' +#' This is a convenience function to set up the lists used as arguments for +#' \code{\link{mkinmod}}. +#' +#' @param submodel Character vector of length one to specify the submodel type. +#' See \code{\link{mkinmod}} for the list of allowed submodel names. +#' @param to Vector of the names of the state variable to which a +#' transformation shall be included in the model. +#' @param sink Should a pathway to sink be included in the model in addition to +#' the pathways to other state variables? +#' @param full_name An optional name to be used e.g. for plotting fits +#' performed with the model. You can use non-ASCII characters here, but then +#' your R code will not be portable, \emph{i.e.} may produce unintended plot +#' results on other operating systems or system configurations. +#' @return A list for use with \code{\link{mkinmod}}. +#' @author Johannes Ranke +#' @examples +#' +#' # One parent compound, one metabolite, both single first order. +#' SFO_SFO <- mkinmod( +#' parent = list(type = "SFO", to = "m1"), +#' m1 = list(type = "SFO")) +#' +#' # The same model using mkinsub +#' SFO_SFO.2 <- mkinmod( +#' parent = mkinsub("SFO", "m1"), +#' m1 = mkinsub("SFO")) +#' +#' # Now supplying full names +#' SFO_SFO.2 <- mkinmod( +#' parent = mkinsub("SFO", "m1", full_name = "Test compound"), +#' m1 = mkinsub("SFO", full_name = "Metabolite M1")) +#' +#' @export mkinsub <- function(submodel, to = NULL, sink = TRUE, full_name = NA) { return(list(type = submodel, to = to, sink = sink, full_name = full_name)) } -# vim: set ts=2 sw=2 expandtab: @@ -1,24 +1,65 @@ -# Copyright (C) 2015,2019 Johannes Ranke -# Contact: jranke@uni-bremen.de -# The summary function is an adapted and extended version of summary.modFit -# from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn -# inspired by summary.nls.lm - -# 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/> - +#' Fit one or more kinetic models with one or more state variables to one or +#' more datasets +#' +#' This function calls \code{\link{mkinfit}} on all combinations of models and +#' datasets specified in its first two arguments. +#' +#' @param models Either a character vector of shorthand names like +#' \code{c("SFO", "FOMC", "DFOP", "HS", "SFORB")}, or an optionally named +#' list of \code{\link{mkinmod}} objects. +#' @param datasets An optionally named list of datasets suitable as observed +#' data for \code{\link{mkinfit}}. +#' @param cores The number of cores to be used for multicore processing. This +#' is only used when the \code{cluster} argument is \code{NULL}. On Windows +#' machines, cores > 1 is not supported, you need to use the \code{cluster} +#' argument to use multiple logical processors. +#' @param cluster A cluster as returned by \code{\link{makeCluster}} to be used +#' for parallel execution. +#' @param \dots Further arguments that will be passed to \code{\link{mkinfit}}. +#' @importFrom parallel mclapply parLapply detectCores +#' @return A matrix of \code{\link{mkinfit}} objects that can be indexed using +#' the model and dataset names as row and column indices. +#' @author Johannes Ranke +#' @seealso \code{\link{[.mmkin}} for subsetting, \code{\link{plot.mmkin}} for +#' plotting. +#' @keywords optimize +#' @examples +#' +#' \dontrun{ +#' m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"), +#' M1 = mkinsub("SFO", "M2"), +#' M2 = mkinsub("SFO"), use_of_ff = "max") +#' +#' m_synth_FOMC_lin <- mkinmod(parent = mkinsub("FOMC", "M1"), +#' M1 = mkinsub("SFO", "M2"), +#' M2 = mkinsub("SFO"), use_of_ff = "max") +#' +#' models <- list(SFO_lin = m_synth_SFO_lin, FOMC_lin = m_synth_FOMC_lin) +#' datasets <- lapply(synthetic_data_for_UBA_2014[1:3], function(x) x$data) +#' names(datasets) <- paste("Dataset", 1:3) +#' +#' time_default <- system.time(fits.0 <- mmkin(models, datasets, quiet = TRUE)) +#' time_1 <- system.time(fits.4 <- mmkin(models, datasets, cores = 1, quiet = TRUE)) +#' +#' time_default +#' time_1 +#' +#' endpoints(fits.0[["SFO_lin", 2]]) +#' +#' # plot.mkinfit handles rows or columns of mmkin result objects +#' plot(fits.0[1, ]) +#' plot(fits.0[1, ], obs_var = c("M1", "M2")) +#' plot(fits.0[, 1]) +#' # Use double brackets to extract a single mkinfit object, which will be plotted +#' # by plot.mkinfit and can be plotted using plot_sep +#' plot(fits.0[[1, 1]], sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE) +#' plot_sep(fits.0[[1, 1]]) +#' # Plotting with mmkin (single brackets, extracting an mmkin object) does not +#' # allow to plot the observed variables separately +#' plot(fits.0[1, 1]) +#' } +#' +#' @export mmkin mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, cores = round(detectCores()/2), cluster = NULL, ...) { @@ -66,7 +107,35 @@ mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, return(results) } -"[.mmkin" <- function(x, i, j, ..., drop = FALSE) { +#' Subsetting method for mmkin objects +#' +#' Subsetting method for mmkin objects. +#' +#' @param x An \code{\link{mmkin} object} +#' @param i Row index selecting the fits for specific models +#' @param j Column index selecting the fits to specific datasets +#' @param ... Not used, only there to satisfy the generic method definition +#' @param drop If FALSE, the method always returns an mmkin object, otherwise +#' either a list of mkinfit objects or a single mkinfit object. +#' @return An object of class \code{\link{mmkin}}. +#' @author Johannes Ranke +#' @rdname Extract.mmkin +#' @examples +#' +#' # Only use one core, to pass R CMD check --as-cran +#' fits <- mmkin(c("SFO", "FOMC"), list(B = FOCUS_2006_B, C = FOCUS_2006_C), +#' cores = 1, quiet = TRUE) +#' fits["FOMC", ] +#' fits[, "B"] +#' fits["SFO", "B"] +#' +#' head( +#' # This extracts an mkinfit object with lots of components +#' fits[["FOMC", "B"]] +#' ) +#' +#' @export +`[.mmkin` <- function(x, i, j, ..., drop = FALSE) { class(x) <- NULL x_sub <- x[i, j, drop = drop] if (!drop) class(x_sub) <- "mmkin" @@ -1,21 +1,40 @@ -# Copyright (C) 2019 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/> - +#' Evaluate parent kinetics using the NAFTA guidance +#' +#' The function fits the SFO, IORE and DFOP models using \code{\link{mmkin}} +#' and returns an object of class \code{nafta} that has methods for printing +#' and plotting. +#' +#' @param ds A dataframe that must contain one variable called "time" with the +#' time values specified by the \code{time} argument, one column called +#' "name" with the grouping of the observed values, and finally one column of +#' observed values called "value". +#' @param title Optional title of the dataset +#' @param quiet Should the evaluation text be shown? +#' @param \dots Further arguments passed to \code{\link{mmkin}} (not for the +#' printing method). +#' @importFrom stats qf +#' @return An list of class \code{nafta}. The list element named "mmkin" is the +#' \code{\link{mmkin}} object containing the fits of the three models. The +#' list element named "title" contains the title of the dataset used. The +#' list element "data" contains the dataset used in the fits. +#' @author Johannes Ranke +#' @source NAFTA (2011) Guidance for evaluating and calculating degradation +#' kinetics in environmental media. NAFTA Technical Working Group on +#' Pesticides +#' \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/guidance-evaluating-and-calculating-degradation} +#' accessed 2019-02-22 +#' +#' US EPA (2015) Standard Operating Procedure for Using the NAFTA Guidance to +#' Calculate Representative Half-life Values and Characterizing Pesticide +#' Degradation +#' \url{https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/standard-operating-procedure-using-nafta-guidance} +#' @examples +#' +#' nafta_evaluation <- nafta(NAFTA_SOP_Appendix_D, cores = 1) +#' print(nafta_evaluation) +#' plot(nafta_evaluation) +#' +#' @export nafta <- function(ds, title = NA, quiet = FALSE, ...) { if (length(levels(ds$name)) > 1) { stop("The NAFTA procedure is only defined for decline data for a single compound") @@ -56,6 +75,20 @@ nafta <- function(ds, title = NA, quiet = FALSE, ...) { return(result) } +#' Plot the results of the three models used in the NAFTA scheme. +#' +#' The plots are ordered with increasing complexity of the model in this +#' function (SFO, then IORE, then DFOP). +#' +#' Calls \code{\link{plot.mmkin}}. +#' +#' @param x An object of class \code{\link{nafta}}. +#' @param legend Should a legend be added? +#' @param main Possibility to override the main title of the plot. +#' @param \dots Further arguments passed to \code{\link{plot.mmkin}}. +#' @return The function is called for its side effect. +#' @author Johannes Ranke +#' @export plot.nafta <- function(x, legend = FALSE, main = "auto", ...) { if (main == "auto") { if (is.na(x$title)) main = "" @@ -64,6 +97,16 @@ plot.nafta <- function(x, legend = FALSE, main = "auto", ...) { plot(x$mmkin, ..., legend = legend, main = main) } +#' Print nafta objects +#' +#' Print nafta objects. The results for the three models are printed in the +#' order of increasing model complexity, i.e. SFO, then IORE, and finally DFOP. +#' +#' @param x An \code{\link{nafta}} object. +#' @param digits Number of digits to be used for printing parameters and +#' dissipation times. +#' @rdname nafta +#' @export print.nafta <- function(x, quiet = TRUE, digits = 3, ...) { cat("Sums of squares:\n") print(x$S) diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index e39da416..16415a3b 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -1,22 +1,88 @@ -# Copyright (C) 2010-2016,2019 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/> if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "observed")) +#' Plot the observed data and the fitted model of an mkinfit object +#' +#' Solves the differential equations with the optimised and fixed parameters +#' from a previous successful call to \code{\link{mkinfit}} and plots the +#' observed data together with the solution of the fitted model. +#' +#' If the current plot device is a \code{\link[tikzDevice]{tikz}} device, then +#' latex is being used for the formatting of the chi2 error level, if +#' \code{show_errmin = TRUE}. +#' +#' @aliases plot.mkinfit plot_sep plot_res plot_err +#' @param x Alias for fit introduced for compatibility with the generic S3 +#' method. +#' @param fit An object of class \code{\link{mkinfit}}. +#' @param obs_vars A character vector of names of the observed variables for +#' which the data and the model should be plotted. Defauls to all observed +#' variables in the model. +#' @param xlab Label for the x axis. +#' @param ylab Label for the y axis. +#' @param xlim Plot range in x direction. +#' @param ylim Plot range in y direction. +#' @param col_obs Colors used for plotting the observed data and the +#' corresponding model prediction lines. +#' @param pch_obs Symbols to be used for plotting the data. +#' @param lty_obs Line types to be used for the model predictions. +#' @param add Should the plot be added to an existing plot? +#' @param legend Should a legend be added to the plot? +#' @param show_residuals Should residuals be shown? If only one plot of the +#' fits is shown, the residual plot is in the lower third of the plot. +#' Otherwise, i.e. if "sep_obs" is given, the residual plots will be located +#' to the right of the plots of the fitted curves. +#' @param show_errplot Should squared residuals and the error model be shown? +#' If only one plot of the fits is shown, this plot is in the lower third of +#' the plot. Otherwise, i.e. if "sep_obs" is given, the residual plots will +#' be located to the right of the plots of the fitted curves. +#' @param maxabs Maximum absolute value of the residuals. This is used for the +#' scaling of the y axis and defaults to "auto". +#' @param sep_obs Should the observed variables be shown in separate subplots? +#' If yes, residual plots requested by "show_residuals" will be shown next +#' to, not below the plot of the fits. +#' @param rel.height.middle The relative height of the middle plot, if more +#' than two rows of plots are shown. +#' @param row_layout Should we use a row layout where the residual plot or the +#' error model plot is shown to the right? +#' @param lpos Position(s) of the legend(s). Passed to \code{\link{legend}} as +#' the first argument. If not length one, this should be of the same length +#' as the obs_var argument. +#' @param inset Passed to \code{\link{legend}} if applicable. +#' @param show_errmin Should the FOCUS chi2 error value be shown in the upper +#' margin of the plot? +#' @param errmin_digits The number of significant digits for rounding the FOCUS +#' chi2 error percentage. +#' @param frame Should a frame be drawn around the plots? +#' @param \dots Further arguments passed to \code{\link{plot}}. +#' @import graphics +#' @importFrom grDevices dev.cur +#' @return The function is called for its side effect. +#' @author Johannes Ranke +#' @examples +#' +#' # One parent compound, one metabolite, both single first order, path from +#' # parent to sink included +#' \dontrun{ +#' SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1", full = "Parent"), +#' m1 = mkinsub("SFO", full = "Metabolite M1" )) +#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE, error_model = "tc") +#' plot(fit) +#' plot_res(fit) +#' plot_err(fit) +#' +#' # Show the observed variables separately, with residuals +#' plot(fit, sep_obs = TRUE, show_residuals = TRUE, lpos = c("topright", "bottomright"), +#' show_errmin = TRUE) +#' +#' # The same can be obtained with less typing, using the convenience function plot_sep +#' plot_sep(fit, lpos = c("topright", "bottomright")) +#' +#' # Show the observed variables separately, with the error model +#' plot(fit, sep_obs = TRUE, show_errplot = TRUE, lpos = c("topright", "bottomright"), +#' show_errmin = TRUE) +#' } +#' +#' @export plot.mkinfit <- function(x, fit = x, obs_vars = names(fit$mkinmod$map), xlab = "Time", ylab = "Observed", @@ -206,17 +272,39 @@ plot.mkinfit <- function(x, fit = x, } if (do_layout) par(oldpar, no.readonly = TRUE) } -# Convenience function for switching on some features of mkinfit -# that have not been made the default to keep compatibility + +#' @rdname plot.mkinfit +#' @export plot_sep <- function(fit, show_errmin = TRUE, ...) { plot.mkinfit(fit, sep_obs = TRUE, show_residuals = TRUE, show_errmin = show_errmin, ...) } + +#' @rdname plot.mkinfit +#' @export plot_res <- function(fit, sep_obs = FALSE, show_errmin = sep_obs, ...) { plot.mkinfit(fit, sep_obs = sep_obs, show_errmin = show_errmin, show_residuals = TRUE, row_layout = TRUE, ...) } + +#' @rdname plot.mkinfit +#' @export plot_err <- function(fit, sep_obs = FALSE, show_errmin = sep_obs, ...) { plot.mkinfit(fit, sep_obs = sep_obs, show_errmin = show_errmin, show_errplot = TRUE, row_layout = TRUE, ...) } + +#' Plot the observed data and the fitted model of an mkinfit object +#' +#' Deprecated function. It now only calls the plot method +#' \code{\link{plot.mkinfit}}. +#' +#' @param fit an object of class \code{\link{mkinfit}}. +#' @param \dots further arguments passed to \code{\link{plot.mkinfit}}. +#' @return The function is called for its side effect. +#' @author Johannes Ranke +#' @export +mkinplot <- function(fit, ...) +{ + plot(fit, ...) +} diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R index ef80949c..eefafe12 100644 --- a/R/plot.mmkin.R +++ b/R/plot.mmkin.R @@ -1,21 +1,53 @@ -# Copyright (C) 2015-2016,2019 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/> - +#' Plot model fits (observed and fitted) and the residuals for a row or column +#' of an mmkin object +#' +#' When x is a row selected from an mmkin object (\code{\link{[.mmkin}}), the +#' same model fitted for at least one dataset is shown. When it is a column, +#' the fit of at least one model to the same dataset is shown. +#' +#' If the current plot device is a \code{\link[tikzDevice]{tikz}} device, then +#' latex is being used for the formatting of the chi2 error level. +#' +#' @param x An object of class \code{\link{mmkin}}, with either one row or one +#' column. +#' @param main The main title placed on the outer margin of the plot. +#' @param legends An index for the fits for which legends should be shown. +#' @param resplot Should the residuals plotted against time, using +#' \code{\link{mkinresplot}}, or as squared residuals against predicted +#' values, with the error model, using \code{\link{mkinerrplot}}. +#' @param show_errmin Should the chi2 error level be shown on top of the plots +#' to the left? +#' @param errmin_var The variable for which the FOCUS chi2 error value should +#' be shown. +#' @param errmin_digits The number of significant digits for rounding the FOCUS +#' chi2 error percentage. +#' @param cex Passed to the plot functions and \code{\link{mtext}}. +#' @param rel.height.middle The relative height of the middle plot, if more +#' than two rows of plots are shown. +#' @param \dots Further arguments passed to \code{\link{plot.mkinfit}} and +#' \code{\link{mkinresplot}}. +#' @return The function is called for its side effect. +#' @author Johannes Ranke +#' @examples +#' +#' \dontrun{ +#' # Only use one core not to offend CRAN checks +#' fits <- mmkin(c("FOMC", "HS"), +#' list("FOCUS B" = FOCUS_2006_B, "FOCUS C" = FOCUS_2006_C), # named list for titles +#' cores = 1, quiet = TRUE, error_model = "tc") +#' plot(fits[, "FOCUS C"]) +#' plot(fits["FOMC", ]) +#' +#' # We can also plot a single fit, if we like the way plot.mmkin works, but then the plot +#' # height should be smaller than the plot width (this is not possible for the html pages +#' # generated by pkgdown, as far as I know). +#' plot(fits["FOMC", "FOCUS C"]) # same as plot(fits[1, 2]) +#' +#' # Show the error models +#' plot(fits["FOMC", ], resplot = "errmod") +#' } +#' +#' @export plot.mmkin <- function(x, main = "auto", legends = 1, resplot = c("time", "errmod"), show_errmin = TRUE, diff --git a/R/sigma_twocomp.R b/R/sigma_twocomp.R index b06816c1..c9a15aa8 100644 --- a/R/sigma_twocomp.R +++ b/R/sigma_twocomp.R @@ -1,3 +1,29 @@ +#' Two component error model +#' +#' Function describing the standard deviation of the measurement error in +#' dependence of the measured value \eqn{y}: +#' +#' \deqn{\sigma = \sqrt{ \sigma_{low}^2 + y^2 * {rsd}_{high}^2}} sigma = +#' sqrt(sigma_low^2 + y^2 * rsd_high^2) +#' +#' This is the error model used for example by Werner et al. (1978). The model +#' proposed by Rocke and Lorenzato (1995) can be written in this form as well, +#' but assumes approximate lognormal distribution of errors for high values of +#' y. +#' +#' @param y The magnitude of the observed value +#' @param sigma_low The asymptotic minimum of the standard deviation for low +#' observed values +#' @param rsd_high The coefficient describing the increase of the standard +#' deviation with the magnitude of the observed value +#' @return The standard deviation of the response variable. +#' @references Werner, Mario, Brooks, Samuel H., and Knott, Lancaster B. (1978) +#' Additive, Multiplicative, and Mixed Analytical Errors. Clinical Chemistry +#' 24(11), 1895-1898. +#' +#' Rocke, David M. and Lorenzato, Stefan (1995) A two-component model for +#' measurement error in analytical chemistry. Technometrics 37(2), 176-184. +#' @export sigma_twocomp <- function(y, sigma_low, rsd_high) { sqrt(sigma_low^2 + y^2 * rsd_high^2) } diff --git a/R/summary.mkinfit.R b/R/summary.mkinfit.R new file mode 100644 index 00000000..90f32da9 --- /dev/null +++ b/R/summary.mkinfit.R @@ -0,0 +1,272 @@ +#' Summary method for class "mkinfit" +#' +#' Lists model equations, initial parameter values, optimised parameters with +#' some uncertainty statistics, the chi2 error levels calculated according to +#' FOCUS guidance (2006) as defined therein, formation fractions, DT50 values +#' and optionally the data, consisting of observed, predicted and residual +#' values. +#' +#' @param object an object of class \code{\link{mkinfit}}. +#' @param x an object of class \code{summary.mkinfit}. +#' @param data logical, indicating whether the data should be included in the +#' summary. +#' @param distimes logical, indicating whether DT50 and DT90 values should be +#' included. +#' @param alpha error level for confidence interval estimation from t +#' distribution +#' @param digits Number of digits to use for printing +#' @param \dots optional arguments passed to methods like \code{print}. +#' @importFrom stats qt pt cov2cor +#' @return The summary function returns a list with components, among others +#' \item{version, Rversion}{The mkin and R versions used} +#' \item{date.fit, date.summary}{The dates where the fit and the summary were +#' produced} +#' \item{diffs}{The differential equations used in the model} +#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions} +#' \item{bpar}{Optimised and backtransformed +#' parameters} +#' \item{data}{The data (see Description above).} +#' \item{start}{The starting values and bounds, if applicable, for optimised +#' parameters.} +#' \item{fixed}{The values of fixed parameters.} +#' \item{errmin }{The chi2 error levels for +#' each observed variable.} +#' \item{bparms.ode}{All backtransformed ODE +#' parameters, for use as starting parameters for related models.} +#' \item{errparms}{Error model parameters.} +#' \item{ff}{The estimated formation fractions derived from the fitted +#' model.} +#' \item{distimes}{The DT50 and DT90 values for each observed variable.} +#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} +#' The print method is called for its side effect, i.e. printing the summary. +#' @author Johannes Ranke +#' @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} +#' @examples +#' +#' summary(mkinfit(mkinmod(parent = mkinsub("SFO")), FOCUS_2006_A, quiet = TRUE)) +#' +#' @export +summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) { + param <- object$par + pnames <- names(param) + bpnames <- names(object$bparms.optim) + epnames <- names(object$errparms) + p <- length(param) + mod_vars <- names(object$mkinmod$diffs) + covar <- try(solve(object$hessian), silent = TRUE) + covar_notrans <- try(solve(object$hessian_notrans), silent = TRUE) + rdf <- object$df.residual + + if (!is.numeric(covar) | is.na(covar[1])) { + covar <- NULL + se <- lci <- uci <- rep(NA, p) + } else { + rownames(covar) <- colnames(covar) <- pnames + se <- sqrt(diag(covar)) + lci <- param + qt(alpha/2, rdf) * se + uci <- param + qt(1-alpha/2, rdf) * se + } + + beparms.optim <- c(object$bparms.optim, object$par[epnames]) + if (!is.numeric(covar_notrans) | is.na(covar_notrans[1])) { + covar_notrans <- NULL + se_notrans <- tval <- pval <- rep(NA, p) + } else { + rownames(covar_notrans) <- colnames(covar_notrans) <- c(bpnames, epnames) + se_notrans <- sqrt(diag(covar_notrans)) + tval <- beparms.optim / se_notrans + pval <- pt(abs(tval), rdf, lower.tail = FALSE) + } + + names(se) <- pnames + + param <- cbind(param, se, lci, uci) + dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper")) + + bparam <- cbind(Estimate = beparms.optim, se_notrans, + "t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA) + + # Transform boundaries of CI for one parameter at a time, + # with the exception of sets of formation fractions (single fractions are OK). + f_names_skip <- character(0) + for (box in mod_vars) { # Figure out sets of fractions to skip + f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE) + n_paths <- length(f_names) + if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names) + } + + for (pname in pnames) { + if (!pname %in% f_names_skip) { + par.lower <- param[pname, "Lower"] + par.upper <- param[pname, "Upper"] + names(par.lower) <- names(par.upper) <- pname + bpl <- backtransform_odeparms(par.lower, object$mkinmod, + object$transform_rates, + object$transform_fractions) + bpu <- backtransform_odeparms(par.upper, object$mkinmod, + object$transform_rates, + object$transform_fractions) + bparam[names(bpl), "Lower"] <- bpl + bparam[names(bpu), "Upper"] <- bpu + } + } + bparam[epnames, c("Lower", "Upper")] <- param[epnames, c("Lower", "Upper")] + + ans <- list( + version = as.character(utils::packageVersion("mkin")), + Rversion = paste(R.version$major, R.version$minor, sep="."), + date.fit = object$date, + date.summary = date(), + solution_type = object$solution_type, + warning = object$warning, + use_of_ff = object$mkinmod$use_of_ff, + error_model_algorithm = object$error_model_algorithm, + df = c(p, rdf), + covar = covar, + covar_notrans = covar_notrans, + err_mod = object$err_mod, + niter = object$iterations, + calls = object$calls, + time = object$time, + par = param, + bpar = bparam) + + if (!is.null(object$version)) { + ans$fit_version <- object$version + ans$fit_Rversion <- object$Rversion + } + + ans$diffs <- object$mkinmod$diffs + if(data) ans$data <- object$data + ans$start <- object$start + ans$start_transformed <- object$start_transformed + + ans$fixed <- object$fixed + + ans$errmin <- mkinerrmin(object, alpha = 0.05) + + if (object$calls > 0) { + if (!is.null(ans$covar)){ + Corr <- cov2cor(ans$covar) + rownames(Corr) <- colnames(Corr) <- rownames(ans$par) + ans$Corr <- Corr + } else { + warning("Could not calculate correlation; no covariance matrix") + } + } + + ans$bparms.ode <- object$bparms.ode + ep <- endpoints(object) + if (length(ep$ff) != 0) + ans$ff <- ep$ff + if (distimes) ans$distimes <- ep$distimes + if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB + if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message + class(ans) <- c("summary.mkinfit", "summary.modFit") + return(ans) +} + +#' @rdname summary.mkinfit +#' @export +print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) { + if (is.null(x$fit_version)) { + cat("mkin version: ", x$version, "\n") + cat("R version: ", x$Rversion, "\n") + } else { + cat("mkin version used for fitting: ", x$fit_version, "\n") + cat("R version used for fitting: ", x$fit_Rversion, "\n") + } + + cat("Date of fit: ", x$date.fit, "\n") + cat("Date of summary:", x$date.summary, "\n") + + if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n") + + cat("\nEquations:\n") + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]]) + writeLines(strwrap(nice_diffs, exdent = 11)) + df <- x$df + rdf <- df[2] + + cat("\nModel predictions using solution type", x$solution_type, "\n") + + cat("\nFitted using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n") + + if (!is.null(x$err_mod)) { + cat("\nError model: ") + cat(switch(x$err_mod, + const = "Constant variance", + obs = "Variance unique to each observed variable", + tc = "Two-component variance function"), "\n") + + cat("\nError model algorithm:", x$error_model_algorithm, "\n") + if (!is.null(x$d_3_message)) cat(x$d_3_message, "\n") + } + + cat("\nStarting values for parameters to be optimised:\n") + print(x$start) + + cat("\nStarting values for the transformed parameters actually optimised:\n") + print(x$start_transformed) + + cat("\nFixed parameter values:\n") + if(length(x$fixed$value) == 0) cat("None\n") + else print(x$fixed) + + cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n") + print(signif(x$par, digits = digits)) + + if (x$calls > 0) { + cat("\nParameter correlation:\n") + if (!is.null(x$covar)){ + print(x$Corr, digits = digits, ...) + } else { + cat("No covariance matrix") + } + } + + cat("\nBacktransformed parameters:\n") + cat("Confidence intervals for internally transformed parameters are asymmetric.\n") + if ((x$version) < "0.9-36") { + cat("To get the usual (questionable) t-test, upgrade mkin and repeat the fit.\n") + print(signif(x$bpar, digits = digits)) + } else { + cat("t-test (unrealistically) based on the assumption of normal distribution\n") + cat("for estimators of untransformed parameters.\n") + print(signif(x$bpar[, c(1, 3, 4, 5, 6)], digits = digits)) + } + + cat("\nFOCUS Chi2 error levels in percent:\n") + x$errmin$err.min <- 100 * x$errmin$err.min + print(x$errmin, digits=digits,...) + + printSFORB <- !is.null(x$SFORB) + if(printSFORB){ + cat("\nEstimated Eigenvalues of SFORB model(s):\n") + print(x$SFORB, digits=digits,...) + } + + printff <- !is.null(x$ff) + if(printff){ + cat("\nResulting formation fractions:\n") + print(data.frame(ff = x$ff), digits=digits,...) + } + + printdistimes <- !is.null(x$distimes) + if(printdistimes){ + cat("\nEstimated disappearance times:\n") + print(x$distimes, digits=digits,...) + } + + printdata <- !is.null(x$data) + if (printdata){ + cat("\nData:\n") + print(format(x$data, digits = digits, ...), row.names = FALSE) + } + + invisible(x) +} diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index f69f4ebd..28e58f87 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -1,180 +1,258 @@ -# Copyright (C) 2010-2014,2019 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/>
-
-transform_odeparms <- function(parms, mkinmod,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # We need the model specification for the names of the model
- # variables and the information on the sink
- spec = mkinmod$spec
-
- # Set up container for transformed parameters
- transparms <- numeric(0)
-
- # Do not transform initial values for state variables
- state.ini.optim <- parms[grep("_0$", names(parms))]
- transparms[names(state.ini.optim)] <- state.ini.optim
-
- # Log transformation for rate constants if requested
- k <- parms[grep("^k_", names(parms))]
- k__iore <- parms[grep("^k__iore_", names(parms))]
- k <- c(k, k__iore)
- if (length(k) > 0) {
- if(transform_rates) {
- transparms[paste0("log_", names(k))] <- log(k)
- } else transparms[names(k)] <- k
- }
-
- # Do not transform exponents in IORE models
- N <- parms[grep("^N", names(parms))]
- transparms[names(N)] <- N
-
- # Go through state variables and apply isometric logratio transformation to
- # formation fractions if requested
- mod_vars = names(spec)
- for (box in mod_vars) {
- f <- parms[grep(paste("^f", box, sep = "_"), names(parms))]
-
- if (length(f) > 0) {
- if(transform_fractions) {
- if (spec[[box]]$sink) {
- trans_f <- ilr(c(f, 1 - sum(f)))
- trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
- transparms[trans_f_names] <- trans_f
- } else {
- if (length(f) > 1) {
- trans_f <- ilr(f)
- trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
- transparms[trans_f_names] <- trans_f
- }
- }
- } else {
- transparms[names(f)] <- f
- }
- }
- }
-
- # Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2
- # and HS parameter tb as well as logistic model parameters kmax, k0 and r if
- # transformation of rates is requested
- for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
- if (!is.na(parms[pname])) {
- if (transform_rates) {
- transparms[paste0("log_", pname)] <- log(parms[pname])
- } else {
- transparms[pname] <- parms[pname]
- }
- }
- }
-
- # DFOP parameter g is treated as a fraction
- if (!is.na(parms["g"])) {
- g <- parms["g"]
- if (transform_fractions) {
- transparms["g_ilr"] <- ilr(c(g, 1 - g))
- } else {
- transparms["g"] <- g
- }
- }
-
- return(transparms)
-}
-
-backtransform_odeparms <- function(transparms, mkinmod,
- transform_rates = TRUE,
- transform_fractions = TRUE)
-{
- # We need the model specification for the names of the model
- # variables and the information on the sink
- spec = mkinmod$spec
-
- # Set up container for backtransformed parameters
- parms <- numeric(0)
-
- # Do not transform initial values for state variables
- state.ini.optim <- transparms[grep("_0$", names(transparms))]
- parms[names(state.ini.optim)] <- state.ini.optim
-
- # Exponential transformation for rate constants
- if(transform_rates) {
- trans_k <- transparms[grep("^log_k_", names(transparms))]
- trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))]
- trans_k = c(trans_k, trans_k__iore)
- if (length(trans_k) > 0) {
- k_names <- gsub("^log_k", "k", names(trans_k))
- parms[k_names] <- exp(trans_k)
- }
- } else {
- trans_k <- transparms[grep("^k_", names(transparms))]
- parms[names(trans_k)] <- trans_k
- trans_k__iore <- transparms[grep("^k__iore_", names(transparms))]
- parms[names(trans_k__iore)] <- trans_k__iore
- }
-
- # Do not transform exponents in IORE models
- N <- transparms[grep("^N", names(transparms))]
- parms[names(N)] <- N
-
- # Go through state variables and apply inverse isometric logratio transformation
- mod_vars = names(spec)
- for (box in mod_vars) {
- # Get the names as used in the model
- f_names = grep(paste("^f", box, sep = "_"), mkinmod$parms, value = TRUE)
- # Get the formation fraction parameters
- trans_f = transparms[grep(paste("^f", box, sep = "_"), names(transparms))]
- if (length(trans_f) > 0) {
- if(transform_fractions) {
- f <- invilr(trans_f)
- if (spec[[box]]$sink) {
- parms[f_names] <- f[1:length(f)-1]
- } else {
- parms[f_names] <- f
- }
- } else {
- parms[names(trans_f)] <- trans_f
- }
- }
- }
-
- # Transform parameters also for FOMC, DFOP, HS and logistic models
- for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
- if (transform_rates) {
- pname_trans = paste0("log_", pname)
- if (!is.na(transparms[pname_trans])) {
- parms[pname] <- exp(transparms[pname_trans])
- }
- } else {
- if (!is.na(transparms[pname])) {
- parms[pname] <- transparms[pname]
- }
- }
- }
-
- # DFOP parameter g is treated as a fraction
- if (!is.na(transparms["g_ilr"])) {
- g_ilr <- transparms["g_ilr"]
- parms["g"] <- invilr(g_ilr)[1]
- }
- if (!is.na(transparms["g"])) {
- parms["g"] <- transparms["g"]
- }
-
- return(parms)
-}
-# vim: set ts=2 sw=2 expandtab:
+#' Functions to transform and backtransform kinetic parameters for fitting +#' +#' The transformations are intended to map parameters that should only take on +#' restricted values to the full scale of real numbers. For kinetic rate +#' constants and other paramters that can only take on positive values, a +#' simple log transformation is used. For compositional parameters, such as the +#' formations fractions that should always sum up to 1 and can not be negative, +#' the \code{\link{ilr}} transformation is used. +#' +#' The transformation of sets of formation fractions is fragile, as it supposes +#' the same ordering of the components in forward and backward transformation. +#' This is no problem for the internal use in \code{\link{mkinfit}}. +#' +#' @aliases transform_odeparms backtransform_odeparms +#' @param parms Parameters of kinetic models as used in the differential +#' equations. +#' @param transparms Transformed parameters of kinetic models as used in the +#' fitting procedure. +#' @param mkinmod The kinetic model of class \code{\link{mkinmod}}, containing +#' the names of the model variables that are needed for grouping the +#' formation fractions before \code{\link{ilr}} transformation, the parameter +#' names and the information if the pathway to sink is included in the model. +#' @param transform_rates Boolean specifying if kinetic rate constants should +#' be transformed in the model specification used in the fitting for better +#' compliance with the assumption of normal distribution of the estimator. If +#' TRUE, also alpha and beta parameters of the FOMC model are +#' log-transformed, as well as k1 and k2 rate constants for the DFOP and HS +#' models and the break point tb of the HS model. +#' @param transform_fractions Boolean specifying if formation fractions +#' constants should be transformed in the model specification used in the +#' fitting for better compliance with the assumption of normal distribution +#' of the estimator. The default (TRUE) is to do transformations. The g +#' parameter of the DFOP and HS models are also transformed, as they can also +#' be seen as compositional data. The transformation used for these +#' transformations is the \code{\link{ilr}} transformation. +#' @return A vector of transformed or backtransformed parameters with the same +#' names as the original parameters. +#' @author Johannes Ranke +#' @examples +#' +#' SFO_SFO <- mkinmod( +#' parent = list(type = "SFO", to = "m1", sink = TRUE), +#' m1 = list(type = "SFO")) +#' # Fit the model to the FOCUS example dataset D using defaults +#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) +#' fit.s <- summary(fit) +#' # Transformed and backtransformed parameters +#' print(fit.s$par, 3) +#' print(fit.s$bpar, 3) +#' +#' \dontrun{ +#' # Compare to the version without transforming rate parameters +#' fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE, quiet = TRUE) +#' fit.2.s <- summary(fit.2) +#' print(fit.2.s$par, 3) +#' print(fit.2.s$bpar, 3) +#' } +#' +#' initials <- fit$start$value +#' names(initials) <- rownames(fit$start) +#' transformed <- fit$start_transformed$value +#' names(transformed) <- rownames(fit$start_transformed) +#' transform_odeparms(initials, SFO_SFO) +#' backtransform_odeparms(transformed, SFO_SFO) +#' +#' \dontrun{ +#' # The case of formation fractions +#' SFO_SFO.ff <- mkinmod( +#' parent = list(type = "SFO", to = "m1", sink = TRUE), +#' m1 = list(type = "SFO"), +#' use_of_ff = "max") +#' +#' fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE) +#' fit.ff.s <- summary(fit.ff) +#' print(fit.ff.s$par, 3) +#' print(fit.ff.s$bpar, 3) +#' initials <- c("f_parent_to_m1" = 0.5) +#' transformed <- transform_odeparms(initials, SFO_SFO.ff) +#' backtransform_odeparms(transformed, SFO_SFO.ff) +#' +#' # And without sink +#' SFO_SFO.ff.2 <- mkinmod( +#' parent = list(type = "SFO", to = "m1", sink = FALSE), +#' m1 = list(type = "SFO"), +#' use_of_ff = "max") +#' +#' +#' fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_2006_D, quiet = TRUE) +#' fit.ff.2.s <- summary(fit.ff.2) +#' print(fit.ff.2.s$par, 3) +#' print(fit.ff.2.s$bpar, 3) +#' } +#' +#' @export transform_odeparms +transform_odeparms <- function(parms, mkinmod, + transform_rates = TRUE, + transform_fractions = TRUE) +{ + # We need the model specification for the names of the model + # variables and the information on the sink + spec = mkinmod$spec + + # Set up container for transformed parameters + transparms <- numeric(0) + + # Do not transform initial values for state variables + state.ini.optim <- parms[grep("_0$", names(parms))] + transparms[names(state.ini.optim)] <- state.ini.optim + + # Log transformation for rate constants if requested + k <- parms[grep("^k_", names(parms))] + k__iore <- parms[grep("^k__iore_", names(parms))] + k <- c(k, k__iore) + if (length(k) > 0) { + if(transform_rates) { + transparms[paste0("log_", names(k))] <- log(k) + } else transparms[names(k)] <- k + } + + # Do not transform exponents in IORE models + N <- parms[grep("^N", names(parms))] + transparms[names(N)] <- N + + # Go through state variables and apply isometric logratio transformation to + # formation fractions if requested + mod_vars = names(spec) + for (box in mod_vars) { + f <- parms[grep(paste("^f", box, sep = "_"), names(parms))] + + if (length(f) > 0) { + if(transform_fractions) { + if (spec[[box]]$sink) { + trans_f <- ilr(c(f, 1 - sum(f))) + trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_") + transparms[trans_f_names] <- trans_f + } else { + if (length(f) > 1) { + trans_f <- ilr(f) + trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_") + transparms[trans_f_names] <- trans_f + } + } + } else { + transparms[names(f)] <- f + } + } + } + + # Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2 + # and HS parameter tb as well as logistic model parameters kmax, k0 and r if + # transformation of rates is requested + for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) { + if (!is.na(parms[pname])) { + if (transform_rates) { + transparms[paste0("log_", pname)] <- log(parms[pname]) + } else { + transparms[pname] <- parms[pname] + } + } + } + + # DFOP parameter g is treated as a fraction + if (!is.na(parms["g"])) { + g <- parms["g"] + if (transform_fractions) { + transparms["g_ilr"] <- ilr(c(g, 1 - g)) + } else { + transparms["g"] <- g + } + } + + return(transparms) +} + +#' @describeIn transform_odeparms Backtransform the set of transformed parameters +#' @export backtransform_odeparms +backtransform_odeparms <- function(transparms, mkinmod, + transform_rates = TRUE, + transform_fractions = TRUE) +{ + # We need the model specification for the names of the model + # variables and the information on the sink + spec = mkinmod$spec + + # Set up container for backtransformed parameters + parms <- numeric(0) + + # Do not transform initial values for state variables + state.ini.optim <- transparms[grep("_0$", names(transparms))] + parms[names(state.ini.optim)] <- state.ini.optim + + # Exponential transformation for rate constants + if(transform_rates) { + trans_k <- transparms[grep("^log_k_", names(transparms))] + trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))] + trans_k = c(trans_k, trans_k__iore) + if (length(trans_k) > 0) { + k_names <- gsub("^log_k", "k", names(trans_k)) + parms[k_names] <- exp(trans_k) + } + } else { + trans_k <- transparms[grep("^k_", names(transparms))] + parms[names(trans_k)] <- trans_k + trans_k__iore <- transparms[grep("^k__iore_", names(transparms))] + parms[names(trans_k__iore)] <- trans_k__iore + } + + # Do not transform exponents in IORE models + N <- transparms[grep("^N", names(transparms))] + parms[names(N)] <- N + + # Go through state variables and apply inverse isometric logratio transformation + mod_vars = names(spec) + for (box in mod_vars) { + # Get the names as used in the model + f_names = grep(paste("^f", box, sep = "_"), mkinmod$parms, value = TRUE) + # Get the formation fraction parameters + trans_f = transparms[grep(paste("^f", box, sep = "_"), names(transparms))] + if (length(trans_f) > 0) { + if(transform_fractions) { + f <- invilr(trans_f) + if (spec[[box]]$sink) { + parms[f_names] <- f[1:length(f)-1] + } else { + parms[f_names] <- f + } + } else { + parms[names(trans_f)] <- trans_f + } + } + } + + # Transform parameters also for FOMC, DFOP, HS and logistic models + for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) { + if (transform_rates) { + pname_trans = paste0("log_", pname) + if (!is.na(transparms[pname_trans])) { + parms[pname] <- exp(transparms[pname_trans]) + } + } else { + if (!is.na(transparms[pname])) { + parms[pname] <- transparms[pname] + } + } + } + + # DFOP parameter g is treated as a fraction + if (!is.na(transparms["g_ilr"])) { + g_ilr <- transparms["g_ilr"] + parms["g"] <- invilr(g_ilr)[1] + } + if (!is.na(transparms["g"])) { + parms["g"] <- transparms["g"] + } + + return(parms) +} +# vim: set ts=2 sw=2 expandtab: |