diff options
37 files changed, 292 insertions, 341 deletions
diff --git a/.Rbuildignore b/.Rbuildignore index 79f6ea09..5ac1f320 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -2,6 +2,7 @@ ^build.log$ ^check.log$ ^test.log$ +^test.R$ ^README.html$ ^mkin.Rcheck ^inst/web/* @@ -20,3 +20,4 @@ vignettes/*_cache/ vignettes/*_files/ mkin.Rcheck install.log +test.R diff --git a/R/SFORB.solution.R b/R/SFORB.solution.R index 45a4533a..4cb94def 100644 --- a/R/SFORB.solution.R +++ b/R/SFORB.solution.R @@ -2,8 +2,8 @@ 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 *
+
+ 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 0995e634..4d998e94 100644 --- a/R/add_err.R +++ b/R/add_err.R @@ -16,8 +16,8 @@ # 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, - n = 1000, LOD = 0.1, reps = 2, +add_err = function(prediction, sdfunc, + n = 1000, LOD = 0.1, reps = 2, digits = 1, seed = NA) { if (!is.na(seed)) set.seed(seed) @@ -32,7 +32,7 @@ add_err = function(prediction, sdfunc, for (i in 1:n) { d_rep = data.frame(lapply(d_long, rep, each = 2)) d_rep$value = rnorm(length(d_rep$value), d_rep$value, sdfunc(d_rep$value)) - + d_rep[d_rep$time == 0 & d_rep$name %in% c("M1", "M2"), "value"] <- 0 # Set values below the LOD to NA diff --git a/R/endpoints.R b/R/endpoints.R index 03c38ee0..ac1e3e7c 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -8,19 +8,19 @@ endpoints <- function(fit) { 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)),
+ 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]
+ 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 = "_"),
+ 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]]
@@ -91,7 +91,7 @@ endpoints <- function(fit) { 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
}
@@ -1,5 +1,3 @@ -# $Id$ - # Copyright (C) 2012 René Lehmann, Johannes Ranke # Contact: jranke@uni-bremen.de @@ -44,7 +42,7 @@ invilr<-function(x) { # Work around a numerical problem with NaN values returned by the above # Only works if there is only one NaN value: replace it with 1 # if the sum of the other components is < 1e-10 - if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10) + if (sum(is.na(z)) == 1 && sum(z[!is.na(z)]) < 1e-10) z = ifelse(is.na(z), 1, z) return(z) diff --git a/R/mkin_long_to_wide.R b/R/mkin_long_to_wide.R index 87e1afd5..081262f8 100644 --- a/R/mkin_long_to_wide.R +++ b/R/mkin_long_to_wide.R @@ -1,7 +1,5 @@ -# $Id$
-
# Copyright (C) 2010-2011 Johannes Ranke
-# Contact: mkin-devel@lists.berlios.de
+# Contact: jranke@uni-bremen.de
# This file is part of the R package mkin
@@ -31,7 +31,7 @@ #' @field data A dataframe with at least the columns name, time and value #' in order to be compatible with mkinfit #' @example inst/examples/mkinds.R -mkinds <- R6Class("mkinds", +mkinds <- R6Class("mkinds", public = list( title = NULL, sampling_times = NULL, diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 6103dd1d..b361c466 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -24,7 +24,7 @@ mkinerrmin <- function(fit, alpha = 0.05) kinerrmin <- function(errdata, n.parms) { means.mean <- mean(errdata$value_mean, na.rm = TRUE) df = length(errdata$value_mean) - n.parms - + err.min <- sqrt((1 / qchisq(1 - alpha, df)) * sum((errdata$value_mean - errdata$value_pred)^2)/(means.mean^2)) @@ -32,12 +32,12 @@ mkinerrmin <- function(fit, alpha = 0.05) } means <- aggregate(value ~ time + name, data = fit$observed, mean, na.rm=TRUE) - errdata <- merge(means, fit$predicted, by = c("time", "name"), + errdata <- merge(means, fit$predicted, by = c("time", "name"), suffixes = c("_mean", "_pred")) errdata <- errdata[order(errdata$time, errdata$name), ] # Remove values at time zero for variables whose value for state.ini is fixed, - # as these will not have any effect in the optimization and should therefore not + # as these will not have any effect in the optimization and should therefore not # be counted as degrees of freedom. fixed_initials = gsub("_0$", "", rownames(subset(fit$fixed, type = "state"))) errdata <- subset(errdata, !(time == 0 & name %in% fixed_initials)) @@ -59,10 +59,10 @@ mkinerrmin <- function(fit, alpha = 0.05) # Rate constants and IORE exponents are attributed to the source variable n.k.optim <- length(grep(paste("^k", obs_var, sep="_"), names(parms.optim))) - n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), + n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"), names(parms.optim))) n.k__iore.optim <- length(grep(paste("^k__iore", obs_var, sep="_"), names(parms.optim))) - n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var, + n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var, sep = "_"), names(parms.optim))) @@ -87,8 +87,8 @@ mkinerrmin <- function(fit, alpha = 0.05) # FOMC, DFOP and HS parameters are only counted if we are looking at the # first variable in the model which is always the source variable if (obs_var == fit$obs_vars[[1]]) { - special_parms = c("alpha", "log_alpha", "beta", "log_beta", - "k1", "log_k1", "k2", "log_k2", + special_parms = c("alpha", "log_alpha", "beta", "log_beta", + "k1", "log_k1", "k2", "log_k2", "g", "g_ilr", "tb", "log_tb") n.optim <- n.optim + length(intersect(special_parms, names(parms.optim))) } diff --git a/R/mkinfit.R b/R/mkinfit.R index 4c4a9717..05986ff9 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -23,7 +23,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) mkinfit <- function(mkinmod, observed,
parms.ini = "auto",
- state.ini = "auto",
+ state.ini = "auto",
fixed_parms = NULL,
fixed_initials = names(mkinmod$diffs)[-1],
from_max_mean = FALSE,
@@ -45,7 +45,7 @@ mkinfit <- function(mkinmod, observed, {
# 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")
+ parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE")
if (class(mkinmod) != "mkinmod") {
presumed_parent_name = observed[which.max(observed$value), "name"]
if (mkinmod[[1]] %in% parent_models_available) {
@@ -55,7 +55,7 @@ mkinfit <- function(mkinmod, observed, } else {
stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ",
paste(parent_models_available, collapse = ", "))
- }
+ }
}
# Check optimisation method and set maximum number of iterations if specified
@@ -83,7 +83,7 @@ mkinfit <- function(mkinmod, observed, # 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)
+ 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)
@@ -138,7 +138,7 @@ mkinfit <- function(mkinmod, observed, 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("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
@@ -282,9 +282,9 @@ mkinfit <- function(mkinmod, observed, transform_fractions = transform_fractions)
# Solve the system with current transformed parameter values
- out <- mkinpredict(mkinmod, parms,
- odeini, outtimes,
- solution_type = solution_type,
+ out <- mkinpredict(mkinmod, parms,
+ odeini, outtimes,
+ solution_type = solution_type,
use_compiled = use_compiled,
method.ode = method.ode,
atol = atol, rtol = rtol, ...)
@@ -303,27 +303,27 @@ mkinfit <- function(mkinmod, observed, outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100)
out_plot <- mkinpredict(mkinmod, parms,
- odeini, outtimes_plot,
- solution_type = solution_type,
+ odeini, outtimes_plot,
+ solution_type = solution_type,
use_compiled = use_compiled,
method.ode = method.ode,
atol = atol, rtol = rtol, ...)
- plot(0, type="n",
+ plot(0, type="n",
xlim = range(observed$time), ylim = c(0, max(observed$value, na.rm=TRUE)),
xlab = "Time", ylab = "Observed")
col_obs <- pch_obs <- 1:length(obs_vars)
lty_obs <- rep(1, length(obs_vars))
names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars
for (obs_var in obs_vars) {
- points(subset(observed, name == obs_var, c(time, value)),
+ points(subset(observed, name == obs_var, c(time, value)),
pch = pch_obs[obs_var], col = col_obs[obs_var])
}
matlines(out_plot$time, out_plot[-1], col = col_obs, lty = lty_obs)
- legend("topright", inset=c(0.05, 0.05), legend=obs_vars,
+ legend("topright", inset=c(0.05, 0.05), legend=obs_vars,
col=col_obs, pch=pch_obs, lty=1:length(pch_obs))
}
-
+
assign("cost.old", mC$model, inherits=TRUE)
}
return(mC)
@@ -343,9 +343,9 @@ mkinfit <- function(mkinmod, observed, odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
# Solve the system with current transformed parameter values
- out <- mkinpredict(mkinmod, odeparms,
- odeini, outtimes,
- solution_type = solution_type,
+ out <- mkinpredict(mkinmod, odeparms,
+ odeini, outtimes,
+ solution_type = solution_type,
use_compiled = use_compiled,
method.ode = method.ode,
atol = atol, rtol = rtol, ...)
@@ -381,8 +381,8 @@ mkinfit <- function(mkinmod, observed, # Do the fit and take the time
fit_time <- system.time({
- fit <- modFit(cost, c(state.ini.optim, transparms.optim),
- method = method.modFit, control = control.modFit,
+ fit <- modFit(cost, c(state.ini.optim, transparms.optim),
+ method = method.modFit, control = control.modFit,
lower = lower, upper = upper, ...)
# Reiterate the fit until convergence of the variance components (IRLS)
@@ -423,7 +423,7 @@ mkinfit <- function(mkinmod, observed, # Check for convergence
if (method.modFit == "Marq") {
if (!fit$info %in% c(1, 2, 3)) {
- fit$warning = paste0("Optimisation by method ", method.modFit,
+ fit$warning = paste0("Optimisation by method ", method.modFit,
" did not converge.\n",
"The message returned by nls.lm is:\n",
fit$message)
@@ -435,10 +435,10 @@ mkinfit <- function(mkinmod, observed, }
if (method.modFit %in% c("Port", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {
if (fit$convergence != 0) {
- fit$warning = paste0("Optimisation by method ", method.modFit,
+ fit$warning = paste0("Optimisation by method ", method.modFit,
" did not converge.\n",
"Convergence code is ", fit$convergence,
- ifelse(is.null(fit$message), "",
+ ifelse(is.null(fit$message), "",
paste0("\nMessage is ", fit$message)))
warning(fit$warning)
} else {
@@ -483,7 +483,7 @@ mkinfit <- function(mkinmod, observed, # Estimate the Hessian for the model cost without parameter transformations
# to make it possible to obtain the usual t-test
# Code ported from FME::modFit
- Jac_notrans <- try(gradient(function(p, ...) cost_notrans(p)$residuals$res,
+ Jac_notrans <- try(gradient(function(p, ...) cost_notrans(p)$residuals$res,
bparms.optim, centered = TRUE), silent = TRUE)
if (inherits(Jac_notrans, "try-error")) {
warning("Calculation of the Jacobian failed for the cost function of the untransformed model.\n",
@@ -494,9 +494,9 @@ mkinfit <- function(mkinmod, observed, }
# Collect initial parameter values in three dataframes
- fit$start <- data.frame(value = c(state.ini.optim,
+ fit$start <- data.frame(value = c(state.ini.optim,
parms.optim))
- fit$start$type = c(rep("state", length(state.ini.optim)),
+ fit$start$type = c(rep("state", length(state.ini.optim)),
rep("deparm", length(parms.optim)))
fit$start_transformed = data.frame(
@@ -505,7 +505,7 @@ mkinfit <- function(mkinmod, observed, upper = upper)
fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))
- fit$fixed$type = c(rep("state", length(state.ini.fixed)),
+ fit$fixed$type = c(rep("state", length(state.ini.fixed)),
rep("deparm", length(parms.fixed)))
# Collect observed, predicted, residuals and weighting
@@ -535,18 +535,18 @@ mkinfit <- function(mkinmod, observed, fit$reweight.max.iter <- reweight.max.iter
# Return different sets of backtransformed parameters for summary and plotting
- fit$bparms.optim <- bparms.optim
+ 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.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$date <- date()
- class(fit) <- c("mkinfit", "modFit")
+ class(fit) <- c("mkinfit", "modFit")
return(fit)
}
@@ -602,36 +602,36 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, 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,
+ 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_rates,
object$transform_fractions)
bparam[names(bpl), "Lower"] <- bpl
bparam[names(bpu), "Upper"] <- bpu
- }
+ }
}
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,
- method.modFit = object$method.modFit,
- warning = object$warning,
- use_of_ff = object$mkinmod$use_of_ff,
+ date.fit = object$date,
+ date.summary = date(),
+ solution_type = object$solution_type,
+ method.modFit = object$method.modFit,
+ warning = object$warning,
+ use_of_ff = object$mkinmod$use_of_ff,
weight.ini = object$weight.ini,
reweight.method = object$reweight.method,
residuals = object$residuals,
residualVariance = resvar,
sigma = sqrt(resvar),
modVariance = modVariance,
- df = c(p, rdf),
+ df = c(p, rdf),
cov.unscaled = covar,
cov.scaled = covar * resvar,
- info = object$info,
+ info = object$info,
niter = object$iterations,
calls = object$calls,
time = object$time,
@@ -654,8 +654,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ans$ff <- ep$ff
if(distimes) ans$distimes <- ep$distimes
if(length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
- class(ans) <- c("summary.mkinfit", "summary.modFit")
- return(ans)
+ class(ans) <- c("summary.mkinfit", "summary.modFit")
+ return(ans)
}
# Expanded from print.summary.modFit
@@ -674,7 +674,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("\nModel predictions using solution type", x$solution_type, "\n")
- cat("\nFitted with method", x$method.modFit,
+ cat("\nFitted with method", x$method.modFit,
"using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n")
cat("\nWeighting:", x$weight.ini)
@@ -691,7 +691,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . 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))
diff --git a/R/mkinmod.R b/R/mkinmod.R index 5ec62bbe..7f7f8587 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -44,7 +44,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb # Do not return a coefficient matrix mat when FOMC, IORE, DFOP or HS is used for the parent {{{
if(spec[[1]]$type %in% c("FOMC", "IORE", "DFOP", "HS")) {
- mat = FALSE
+ mat = FALSE
} else mat = TRUE
#}}}
@@ -60,7 +60,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb if(!spec[[varname]]$type %in% c("SFO", "FOMC", "IORE", "DFOP", "HS", "SFORB")) stop(
"Available types are SFO, FOMC, IORE, DFOP, HS and SFORB only")
if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) {
- stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,",
+ stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,",
"which is assumed to be the source compartment"))
}
#}}}
@@ -109,7 +109,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb decline_term <- paste0(decline_term, "^", N)
}
} else { # otherwise no decline term needed here
- decline_term = "0"
+ decline_term = "0"
}
} else {
k_compound <- paste("k", box_1, sep = "_")
@@ -156,7 +156,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } 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
+ # 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 = "_")
@@ -189,9 +189,9 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb {
k_from_to <- paste("k", origin_box, target_box, sep = "_")
parms <- c(parms, k_from_to)
- diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
+ diffs[[origin_box]] <- paste(diffs[[origin_box]], "-",
k_from_to, "*", origin_box)
- diffs[[target_box]] <- paste(diffs[[target_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
@@ -201,7 +201,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb } else {
fraction_to_target = paste("f", origin_box, "to", target, sep = "_")
parms <- c(parms, fraction_to_target)
- diffs[[target_box]] <- paste(diffs[[target_box]], "+",
+ diffs[[target_box]] <- paste(diffs[[target_box]], "+",
fraction_to_target, "*", decline_term)
}
}
@@ -261,7 +261,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb 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 = ""),
+ 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) {
@@ -317,7 +317,7 @@ mkinmod <- function(..., use_of_ff = "min", speclist = NULL, quiet = FALSE, verb 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)
@@ -369,6 +369,6 @@ print.mkinmod <- function(x, ...) { cat("\n")
}
if (is.matrix(x$coefmat)) cat("Coefficient matrix $coefmat available\n")
- if (!is.null(x$cf)) cat("Compiled model $cf available\n")
+ if (!is.null(x$cf)) cat("Compiled model $cf available\n")
}
# vim: set foldmethod=marker ts=2 sw=2 expandtab:
diff --git a/R/mkinparplot.R b/R/mkinparplot.R index 0e8bdf5a..af28e3a8 100644 --- a/R/mkinparplot.R +++ b/R/mkinparplot.R @@ -23,8 +23,8 @@ mkinparplot <- function(object) { if ("g" %in% deparms.optim) fractions.optim <- c("g", fractions.optim) rates.optim.unsorted = setdiff(deparms.optim, union(fractions.optim, N.optim)) rates.optim <- rownames(object$start[rates.optim.unsorted, ]) - n.plot <- c(state.optim = length(state.optim), - rates.optim = length(rates.optim), + n.plot <- c(state.optim = length(state.optim), + rates.optim = length(rates.optim), N.optim = length(N.optim), fractions.optim = length(fractions.optim)) n.plot <- n.plot[n.plot > 0] @@ -41,31 +41,31 @@ mkinparplot <- function(object) { values <- bpar[parnames] values_with_confints <- data.frame(t(subset(data.frame(t(values)), !is.na("Lower")))) xlim = switch(type, - state.optim = range(c(0, unlist(values)), + state.optim = range(c(0, unlist(values)), na.rm = TRUE, finite = TRUE), - rates.optim = range(c(0, unlist(values)), + rates.optim = range(c(0, unlist(values)), na.rm = TRUE, finite = TRUE), - N.optim = range(c(0, 1, unlist(values)), + N.optim = range(c(0, 1, unlist(values)), na.rm = TRUE, finite = TRUE), - fractions.optim = range(c(0, 1, unlist(values)), + fractions.optim = range(c(0, 1, unlist(values)), na.rm = TRUE, finite = TRUE)) parname_index <- length(parnames):1 # Reverse order for strip chart - stripchart(values["Estimate", ][parname_index], + stripchart(values["Estimate", ][parname_index], xlim = xlim, ylim = c(0.5, length(get(type)) + 0.5), yaxt = "n") if (type %in% c("rates.optim", "fractions.optim")) abline(v = 0, lty = 2) if (type %in% c("N.optim", "fractions.optim")) abline(v = 1, lty = 2) position <- ifelse(values["Estimate", ] < mean(xlim), "right", "left") - text(ifelse(position == "left", min(xlim), max(xlim)), - parname_index, parnames, + text(ifelse(position == "left", min(xlim), max(xlim)), + parname_index, parnames, pos = ifelse(position == "left", 4, 2)) values.upper.nonInf <- ifelse(values["Upper", ] == Inf, 1.5 * xlim[[2]], values["Upper", ]) # Suppress warnings for non-existing arrow lengths - suppressWarnings(arrows(as.numeric(values["Lower", ]), parname_index, - as.numeric(values.upper.nonInf), parname_index, + suppressWarnings(arrows(as.numeric(values["Lower", ]), parname_index, + as.numeric(values.upper.nonInf), parname_index, code = 3, angle = 90, length = 0.05)) } par(oldpars) diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 5f11c35a..994af703 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -17,8 +17,8 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see <http://www.gnu.org/licenses/> -mkinpredict <- function(mkinmod, odeparms, odeini, - outtimes, solution_type = "deSolve", +mkinpredict <- function(mkinmod, odeparms, odeini, + outtimes, solution_type = "deSolve", use_compiled = "auto", method.ode = "lsoda", atol = 1e-8, rtol = 1e-10, map_output = TRUE, ...) { @@ -40,12 +40,12 @@ mkinpredict <- function(mkinmod, odeparms, odeini, # Create a function calculating the differentials specified by the model # if necessary if (solution_type == "analytical") { - parent.type = names(mkinmod$map[[1]])[1] + parent.type = names(mkinmod$map[[1]])[1] parent.name = names(mkinmod$diffs)[[1]] o <- switch(parent.type, - SFO = SFO.solution(outtimes, + SFO = SFO.solution(outtimes, evalparse(parent.name), - ifelse(mkinmod$use_of_ff == "min", + ifelse(mkinmod$use_of_ff == "min", evalparse(paste("k", parent.name, "sink", sep="_")), evalparse(paste("k", parent.name, sep="_")))), FOMC = FOMC.solution(outtimes, @@ -53,7 +53,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, evalparse("alpha"), evalparse("beta")), IORE = IORE.solution(outtimes, evalparse(parent.name), - ifelse(mkinmod$use_of_ff == "min", + ifelse(mkinmod$use_of_ff == "min", evalparse(paste("k__iore", parent.name, "sink", sep="_")), evalparse(paste("k__iore", parent.name, sep="_"))), evalparse("N_parent")), @@ -75,18 +75,18 @@ mkinpredict <- function(mkinmod, odeparms, odeini, names(out) <- c("time", sub("_free", "", parent.name)) } if (solution_type == "eigen") { - coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), + coefmat.num <- matrix(sapply(as.vector(mkinmod$coefmat), evalparse), nrow = length(mod_vars)) e <- eigen(coefmat.num) c <- solve(e$vectors, odeini) f.out <- function(t) { e$vectors %*% diag(exp(e$values * t), nrow=length(mod_vars)) %*% c } - o <- matrix(mapply(f.out, outtimes), + o <- matrix(mapply(f.out, outtimes), nrow = length(mod_vars), ncol = length(outtimes)) out <- data.frame(outtimes, t(o)) names(out) <- c("time", mod_vars) - } + } if (solution_type == "deSolve") { if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) { out <- ode( @@ -108,7 +108,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, diffs <- vector() for (box in names(mkinmod$diffs)) { - diffname <- paste("d", box, sep="_") + diffname <- paste("d", box, sep="_") diffs[diffname] <- with(as.list(c(time, state, parms)), eval(parse(text=mkinmod$diffs[[box]]))) } @@ -117,14 +117,14 @@ mkinpredict <- function(mkinmod, odeparms, odeini, out <- ode( y = odeini, times = outtimes, - func = mkindiff, + func = mkindiff, parms = odeparms, method = method.ode, atol = atol, rtol = rtol, ... ) - } + } if (sum(is.na(out)) > 0) { stop("Differential equations were not integrated for all output times because\n", "NaN values occurred in output from ode()") @@ -140,7 +140,7 @@ mkinpredict <- function(mkinmod, odeparms, odeini, out_mapped[var] <- rowSums(out[, mkinmod$map[[var]]]) } } - return(out_mapped) + return(out_mapped) } else { return(out) } diff --git a/R/mkinresplot.R b/R/mkinresplot.R index 82ffd2cb..3650ef4b 100644 --- a/R/mkinresplot.R +++ b/R/mkinresplot.R @@ -17,16 +17,16 @@ # this program. If not, see <http://www.gnu.org/licenses/>
if(getRversion() >= '2.15.1') utils::globalVariables(c("variable", "residual"))
-mkinresplot <- function (object,
+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", ...)
+ maxabs = "auto", legend= TRUE, lpos = "topright", ...)
{
obs_vars_all <- as.character(unique(object$data$variable))
if (length(obs_vars) > 0){
- obs_vars <- intersect(obs_vars_all, obs_vars)
+ obs_vars <- intersect(obs_vars_all, obs_vars)
} else obs_vars <- obs_vars_all
residuals <- subset(object$data, variable %in% obs_vars, residual)
@@ -37,7 +37,7 @@ mkinresplot <- function (object, names(col_obs) <- names(pch_obs) <- obs_vars
plot(0, type = "n",
- xlab = xlab, ylab = ylab,
+ xlab = xlab, ylab = ylab,
xlim = xlim,
ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...)
@@ -49,7 +49,7 @@ mkinresplot <- function (object, abline(h = 0, lty = 2)
if (legend == TRUE) {
- legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
+ legend(lpos, inset = c(0.05, 0.05), legend = obs_vars,
col = col_obs[obs_vars], pch = pch_obs[obs_vars])
}
}
@@ -19,10 +19,10 @@ # You should have received a copy of the GNU General Public License along with # this program. If not, see <http://www.gnu.org/licenses/> -mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, - cores = round(detectCores()/2), cluster = NULL, ...) +mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, + cores = round(detectCores()/2), cluster = NULL, ...) { - parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE") + parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE") n.m <- length(models) n.d <- length(datasets) n.fits <- n.m * n.d @@ -32,10 +32,10 @@ mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, if (!all(sapply(models, function(x) inherits(x, "mkinmod")))) { if (!all(models %in% parent_models_available)) { stop("Please supply models as a list of mkinmod objects or a vector combined of\n ", - paste(parent_models_available, collapse = ", ")) + paste(parent_models_available, collapse = ", ")) } else { names(models) <- models - } + } } else { if (is.null(names(models))) names(models) <- as.character(1:n.m) } diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 00598737..3df6c81c 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -20,15 +20,15 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "variable", "obse plot.mkinfit <- function(x, fit = x, obs_vars = names(fit$mkinmod$map), xlab = "Time", ylab = "Observed", - xlim = range(fit$data$time), + xlim = range(fit$data$time), ylim = "default", col_obs = 1:length(obs_vars), - pch_obs = col_obs, + pch_obs = col_obs, lty_obs = rep(1, length(obs_vars)), - add = FALSE, legend = !add, + add = FALSE, legend = !add, show_residuals = FALSE, maxabs = "auto", sep_obs = FALSE, rel.height.middle = 0.9, - lpos = "topright", inset = c(0.05, 0.05), + lpos = "topright", inset = c(0.05, 0.05), show_errmin = FALSE, errmin_digits = 3, ...) { if (add && show_residuals) stop("If adding to an existing plot we can not show residuals") @@ -53,12 +53,12 @@ plot.mkinfit <- function(x, fit = x, rownames(subset(fit$fixed, type == "deparm"))) odeparms <- parms.all[odenames] - out <- try(mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, + out <- try(mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, solution_type = solution_type, atol = fit$atol, rtol = fit$rtol), silent = TRUE) if (inherits(out, "try-error")) { - out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, + out <- mkinpredict(fit$mkinmod, odeparms, odeini, outtimes, solution_type = solution_type, atol = fit$atol, rtol = fit$rtol, use_compiled = FALSE) } @@ -85,7 +85,7 @@ plot.mkinfit <- function(x, fit = x, # and the middle plots (if n_plot_rows >2) are smaller by rel.height.middle rel.heights <- if (n_plot_rows > 2) c(1, rep(rel.height.middle, n_plot_rows - 2), 1) else rep(1, n_plot_rows) - layout_matrix = matrix(1:n_plots, + layout_matrix = matrix(1:n_plots, n_plot_rows, n_plot_cols, byrow = TRUE) layout(layout_matrix, heights = rel.heights) } else { # else show residuals in the lower third to keep compatibility @@ -104,8 +104,8 @@ plot.mkinfit <- function(x, fit = x, # Set ylim to sensible default, or to the specified value if (ylim[[1]] == "default") { - ylim_row = c(0, max(c(subset(fit$data, variable %in% row_obs_vars)$observed, - subset(fit$data, variable %in% row_obs_vars)$fitted), + ylim_row = c(0, max(c(subset(fit$data, variable %in% row_obs_vars)$observed, + subset(fit$data, variable %in% row_obs_vars)$fitted), na.rm = TRUE)) } else { ylim_row = ylim @@ -115,12 +115,12 @@ plot.mkinfit <- function(x, fit = x, # Margins for top row of plots when we have more than one row # Reduce bottom margin by 2.1 - hides x axis legend if (plot_row == 1 & n_plot_rows > 1) { - par(mar = c(3.0, 4.1, 4.1, 2.1)) + par(mar = c(3.0, 4.1, 4.1, 2.1)) } # Margins for middle rows of plots, if any if (plot_row > 1 & plot_row < n_plot_rows) { - # Reduce top margin by 2 after the first plot as we have no main title, + # Reduce top margin by 2 after the first plot as we have no main title, # reduced plot height, therefore we need rel.height.middle in the layout par(mar = c(3.0, 4.1, 2.1, 2.1)) } @@ -134,14 +134,14 @@ plot.mkinfit <- function(x, fit = x, # Set up the main plot if not to be added to an existing plot if (add == FALSE) { - plot(0, type="n", + plot(0, type="n", xlim = xlim, ylim = ylim_row, xlab = xlab, ylab = ylab, ...) } # Plot the data for (obs_var in row_obs_vars) { - points(subset(fit$data, variable == obs_var, c(time, observed)), + points(subset(fit$data, variable == obs_var, c(time, observed)), pch = pch_obs[obs_var], col = col_obs[obs_var]) } @@ -181,8 +181,8 @@ plot.mkinfit <- function(x, fit = x, residuals <- subset(fit$data, variable %in% row_obs_vars, residual) if (maxabs == "auto") maxabs = max(abs(residuals), na.rm = TRUE) if (!sep_obs) par(mar = c(5, 4, 0, 2) + 0.1) - plot(0, type="n", - xlim = xlim, + plot(0, type="n", + xlim = xlim, ylim = c(-1.2 * maxabs, 1.2 * maxabs), xlab = xlab, ylab = "Residuals") for(obs_var in row_obs_vars){ @@ -197,6 +197,6 @@ plot.mkinfit <- function(x, fit = x, # Convenience function for switching on some features of mkinfit # that have not been made the default to keep compatibility plot_sep <- function(fit, sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE, ...) { - plot.mkinfit(fit, sep_obs = TRUE, show_residuals = TRUE, + plot.mkinfit(fit, sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE, ...) } diff --git a/R/plot.mmkin.R b/R/plot.mmkin.R index 7de91e3e..562bbb71 100644 --- a/R/plot.mmkin.R +++ b/R/plot.mmkin.R @@ -16,12 +16,12 @@ # 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.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", errmin_digits = 3, +plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", errmin_digits = 3, cex = 0.7, rel.height.middle = 0.9, ...) { n.m <- nrow(x) n.d <- ncol(x) - # We can handle either a row (different models, same dataset) + # We can handle either a row (different models, same dataset) # or a column (same model, different datasets) if (n.m > 1 & n.d > 1) stop("Please select fits either for one model or for one dataset") if (n.m == 1 & n.d == 1) loop_over = "none" @@ -53,12 +53,12 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e # Margins for top row of plots when we have more than one row # Reduce bottom margin by 2.1 - hides x axis legend if (i.fit == 1 & n.fits > 1) { - par(mar = c(3.0, 4.1, 4.1, 2.1)) + par(mar = c(3.0, 4.1, 4.1, 2.1)) } # Margins for middle rows of plots, if any if (i.fit > 1 & i.fit < n.fits) { - # Reduce top margin by 2 after the first plot as we have no main title, + # Reduce top margin by 2 after the first plot as we have no main title, # reduced plot height, therefore we need rel.height.middle in the layout par(mar = c(3.0, 4.1, 2.1, 2.1)) } @@ -77,7 +77,7 @@ plot.mmkin <- function(x, main = "auto", legends = 1, errmin_var = "All data", e fit_name <- switch(loop_over, models = rownames(x)[i.fit], datasets = colnames(x)[i.fit], - none = "") + none = "") chi2 <- paste0(signif(100 * mkinerrmin(fit)[errmin_var, "err.min"], errmin_digits), "%") mtext(bquote(.(fit_name) ~ chi^2 ~ "error level" == .(chi2)), cex = cex, line = 0.4) diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 0946ff14..c871c52a 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -17,10 +17,10 @@ # this program. If not, see <http://www.gnu.org/licenses/>
transform_odeparms <- function(parms, mkinmod,
- transform_rates = TRUE,
- transform_fractions = TRUE)
+ transform_rates = TRUE,
+ transform_fractions = TRUE)
{
- # We need the model specification for the names of the model
+ # We need the model specification for the names of the model
# variables and the information on the sink
spec = mkinmod$spec
@@ -79,7 +79,7 @@ transform_odeparms <- function(parms, mkinmod, } else {
transparms[pname] <- parms[pname]
}
- }
+ }
}
# DFOP parameter g is treated as a fraction
@@ -99,7 +99,7 @@ backtransform_odeparms <- function(transparms, mkinmod, transform_rates = TRUE,
transform_fractions = TRUE)
{
- # We need the model specification for the names of the model
+ # We need the model specification for the names of the model
# variables and the information on the sink
spec = mkinmod$spec
@@ -1,6 +1,4 @@ TODO for version 1.0 -- Think about what a user would expect from version 1.0 -- Complete the main package vignette named mkin to include a method description - Improve order of parameters in output Nice to have: diff --git a/man/Extract.mmkin.Rd b/man/Extract.mmkin.Rd index cd40e247..6e0c84e3 100644 --- a/man/Extract.mmkin.Rd +++ b/man/Extract.mmkin.Rd @@ -5,7 +5,7 @@ \method{[}{mmkin}(x, i, j, ..., drop = FALSE) } \description{ - Subsetting method for mmkin objects. + Subsetting method for mmkin objects. } \arguments{ \item{x}{An \code{\link{mmkin} object}} @@ -16,7 +16,7 @@ \item{...}{Not used, only there to satisfy the generic method definition} -\item{drop}{If FALSE, the method always returns an mmkin object, otherwise either +\item{drop}{If FALSE, the method always returns an mmkin object, otherwise either a list of mkinfit objects or a single mkinfit object.} } \value{ @@ -27,7 +27,7 @@ } \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), + fits <- mmkin(c("SFO", "FOMC"), list(B = FOCUS_2006_B, C = FOCUS_2006_C), cores = 1, quiet = TRUE) fits["FOMC", ] fits[, "B"] diff --git a/man/add_err.Rd b/man/add_err.Rd index a0068b59..808edbc6 100644 --- a/man/add_err.Rd +++ b/man/add_err.Rd @@ -10,7 +10,7 @@ } \usage{ add_err(prediction, sdfunc, - n = 1000, LOD = 0.1, reps = 2, + n = 1000, LOD = 0.1, reps = 2, digits = 1, seed = NA) } \arguments{ @@ -36,13 +36,13 @@ The number of digits to which the values will be rounded. } \item{seed}{ - The seed used for the generation of random numbers. If NA, the seed + The seed used for the generation of random numbers. If NA, the seed is not set. } } \value{ A list of datasets compatible with \code{\link{mmkin}}, i.e. - the components of the list are datasets compatible with + the components of the list are datasets compatible with \code{\link{mkinfit}}. } \references{ @@ -80,7 +80,7 @@ names(d_SFO_SFO_err) <- paste("Dataset", 1:3) # Be quiet and use the faster Levenberg-Marquardt algorithm, as the datasets # are easy and examples are run often. Use only one core not to offend CRAN # checks -f_SFO_SFO <- mmkin(list("SFO-SFO" = m_SFO_SFO), +f_SFO_SFO <- mmkin(list("SFO-SFO" = m_SFO_SFO), d_SFO_SFO_err, cores = 1, quiet = TRUE, method.modFit = "Marq") diff --git a/man/endpoints.Rd b/man/endpoints.Rd index ad85622d..08773c3e 100644 --- a/man/endpoints.Rd +++ b/man/endpoints.Rd @@ -5,8 +5,8 @@ Function to calculate endpoints for further use from kinetic models fitted with } \description{ 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 +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. } \usage{ @@ -18,7 +18,7 @@ } } \value{ - The result of the forward or backward transformation. The returned components always + The result of the forward or backward transformation. The returned components always sum to 1 for the case of the inverse log-ratio transformation. } \references{ diff --git a/man/mkin_long_to_wide.Rd b/man/mkin_long_to_wide.Rd index 7f01b7f1..5dd335de 100644 --- a/man/mkin_long_to_wide.Rd +++ b/man/mkin_long_to_wide.Rd @@ -7,14 +7,14 @@ mkin_long_to_wide(long_data, time = "time", outtime = "time") } \description{ - This function takes a dataframe in the long form as required by \code{\link{modCost}} - and converts it into a dataframe with one independent variable and several + This function takes a dataframe in the long form as required by \code{\link{modCost}} + and converts it into a dataframe with one independent variable and several dependent variables as columns. } \arguments{ \item{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 + 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". } \item{time}{ diff --git a/man/mkin_wide_to_long.Rd b/man/mkin_wide_to_long.Rd index cadfe4a2..b531c41c 100644 --- a/man/mkin_wide_to_long.Rd +++ b/man/mkin_wide_to_long.Rd @@ -7,12 +7,12 @@ mkin_wide_to_long(wide_data, time = "t") } \description{ - This function simply takes a dataframe with one independent variable and several + 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{modCost}}. } \arguments{ \item{wide_data}{ - The dataframe must contain one variable with the time values specified by the + 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. } \item{time}{ diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index db2f7fda..f2df212f 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -5,11 +5,11 @@ } \description{ This function uses the Flexible Modelling Environment package - \code{\link{FME}} to create a function calculating the model cost, i.e. the + \code{\link{FME}} to create a function calculating the model cost, i.e. the deviation between the kinetic model and the observed data. This model cost is - then minimised using the Port algorithm \code{\link{nlminb}}, + then minimised using the Port algorithm \code{\link{nlminb}}, using the specified initial or fixed parameters and starting values. - Per default, parameters in the kinetic models are internally transformed in order + Per default, parameters in the kinetic models are internally transformed in order to better satisfy the assumption of a normal distribution of their estimators. In each step of the optimsation, the kinetic model is solved using the function \code{\link{mkinpredict}}. The variance of the residuals for each @@ -17,10 +17,10 @@ using the argument \code{reweight.method = "obs"}. } \usage{ -mkinfit(mkinmod, observed, +mkinfit(mkinmod, observed, parms.ini = "auto", state.ini = "auto", - fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], + fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], from_max_mean = FALSE, solution_type = c("auto", "analytical", "eigen", "deSolve"), method.ode = "lsoda", @@ -30,9 +30,9 @@ mkinfit(mkinmod, observed, control.modFit = list(), transform_rates = TRUE, transform_fractions = TRUE, - plot = FALSE, quiet = FALSE, err = NULL, weight = "none", - scaleVar = FALSE, - atol = 1e-8, rtol = 1e-10, n.outtimes = 100, + plot = FALSE, quiet = FALSE, err = NULL, weight = "none", + scaleVar = FALSE, + atol = 1e-8, rtol = 1e-10, n.outtimes = 100, reweight.method = NULL, reweight.tol = 1e-8, reweight.max.iter = 10, trace_parms = FALSE, ...) @@ -112,12 +112,12 @@ mkinfit(mkinmod, observed, "lsoda" is performant, but sometimes fails to converge. } \item{use_compiled}{ - If set to \code{FALSE}, no compiled version of the \code{\link{mkinmod}} + 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 verion is present. + a compiled verion is present. } \item{method.modFit}{ - The optimisation method passed to \code{\link{modFit}}. + The optimisation method passed to \code{\link{modFit}}. In order to optimally deal with problems where local minima occur, the "Port" algorithm is now used per default as it is less prone to get trapped @@ -138,18 +138,18 @@ mkinfit(mkinmod, observed, } \item{maxit.modFit}{ Maximum number of iterations in the optimisation. If not "auto", this will - be passed to the method called by \code{\link{modFit}}, overriding + be passed to the method called by \code{\link{modFit}}, overriding what may be specified in the next argument \code{control.modFit}. } \item{control.modFit}{ Additional arguments passed to the optimisation method used by - \code{\link{modFit}}. + \code{\link{modFit}}. } \item{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 + 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. @@ -157,7 +157,7 @@ mkinfit(mkinmod, observed, \item{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 + 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 @@ -193,16 +193,16 @@ mkinfit(mkinmod, observed, \item{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 numerical solver if that is used (see \code{solution_type} argument. The default value is 100. } \item{reweight.method}{ The method used for iteratively reweighting residuals, also known as iteratively reweighted least squares (IRLS). Default is NULL, the other method implemented is called "obs", meaning that each - observed variable is assumed to have its own variance, this is + observed variable is assumed to have its own variance, this is estimated from the fit and used for weighting the residuals - in each iteration until convergence of this estimate up to + in each iteration until convergence of this estimate up to \code{reweight.tol} or up to the maximum number of iterations specified by \code{reweight.max.iter}. } @@ -217,19 +217,19 @@ mkinfit(mkinmod, observed, Should a trace of the parameter values be listed? } \item{\dots}{ - Further arguments that will be passed to \code{\link{modFit}}. + Further arguments that will be passed to \code{\link{modFit}}. } } \value{ - A list with "mkinfit" and "modFit" in the class attribute. - A summary can be obtained by \code{\link{summary.mkinfit}}. + A list with "mkinfit" and "modFit" in the class attribute. + A summary can be obtained by \code{\link{summary.mkinfit}}. } \seealso{ Plotting methods \code{\link{plot.mkinfit}} and - \code{\link{mkinparplot}}. + \code{\link{mkinparplot}}. - Fitting of several models to several datasets in a single call to - \code{\link{mmkin}}. + Fitting of several models to several datasets in a single call to + \code{\link{mmkin}}. } \note{ The implementation of iteratively reweighted least squares is inspired by the @@ -238,7 +238,7 @@ mkinfit(mkinmod, observed, other GUI derivative of mkin, sponsored by Syngenta. } \note{ - When using the "IORE" submodel for metabolites, fitting with + 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. @@ -257,13 +257,13 @@ 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, +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, +print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "deSolve"))) coef(fit.deSolve) endpoints(fit.deSolve) @@ -278,7 +278,7 @@ FOMC_SFO <- mkinmod( fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D) # Use starting parameters from parent only FOMC fit fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, plot=TRUE) -fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, +fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, parms.ini = fit.FOMC$bparms.ode, plot=TRUE) # Use stepwise fitting, using optimised parameters from parent only fit, SFORB diff --git a/man/mkinpredict.Rd b/man/mkinpredict.Rd index f9eec532..524abbb5 100644 --- a/man/mkinpredict.Rd +++ b/man/mkinpredict.Rd @@ -43,8 +43,8 @@ "lsoda" is performant, but sometimes fails to converge. } \item{use_compiled}{ - If set to \code{FALSE}, no compiled version of the \code{\link{mkinmod}} - model is used, even if is present. + If set to \code{FALSE}, no compiled version of the \code{\link{mkinmod}} + model is used, even if is present. } \item{atol}{ Absolute error tolerance, passed to \code{\link{ode}}. Default is 1e-8, @@ -56,7 +56,7 @@ } \item{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). + variables (default) or for all state variables (if set to FALSE). } \item{\dots}{ Further arguments passed to the ode solver in case such a solver is used. @@ -71,20 +71,20 @@ \examples{ SFO <- mkinmod(degradinol = list(type = "SFO")) # Compare solution types - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), 0:20, + 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, + 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, + 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, + 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, + 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, + 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,] @@ -94,9 +94,9 @@ # The number of output times used to make a lot of difference until the # default for atol was adjusted - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), + mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), seq(0, 20, by = 0.1))[201,] - mkinpredict(SFO, c(k_degradinol_sink = 0.3), c(degradinol = 100), + mkinpredict(SFO, c(k_degradinol_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! @@ -108,11 +108,11 @@ 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), + 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), + c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), solution_type = "deSolve", use_compiled = FALSE)[201,])) } \keyword{ manip } diff --git a/man/mkinsub.Rd b/man/mkinsub.Rd index 352c4097..84e38e21 100644 --- a/man/mkinsub.Rd +++ b/man/mkinsub.Rd @@ -4,7 +4,7 @@ Function to set up a kinetic submodel for one state variable } \description{ - This is a convenience function to set up the lists used as arguments for + This is a convenience function to set up the lists used as arguments for \code{\link{mkinmod}}. } \usage{ @@ -12,15 +12,15 @@ mkinsub(submodel, to = NULL, sink = TRUE, full_name = NA) } \arguments{ \item{submodel}{ - Character vector of length one to specify the submodel type. See + Character vector of length one to specify the submodel type. See \code{\link{mkinmod}} for the list of allowed submodel names. } \item{to}{ - Vector of the names of the state variable to which a transformation + Vector of the names of the state variable to which a transformation shall be included in the model. } \item{sink}{ - Should a pathway to sink be included in the model in addition to the + Should a pathway to sink be included in the model in addition to the pathways to other state variables? } \item{full_name}{ diff --git a/man/mmkin.Rd b/man/mmkin.Rd index 689a1eb4..3a8d4a1f 100644 --- a/man/mmkin.Rd +++ b/man/mmkin.Rd @@ -26,11 +26,11 @@ mmkin(models, datasets, used when the \code{cluster} argument is \code{NULL}. } \item{cluster}{ - A cluster as returned by \code{\link{makeCluster}} to be used for parallel + A cluster as returned by \code{\link{makeCluster}} to be used for parallel execution. } \item{\dots}{ - Further arguments that will be passed to \code{\link{mkinfit}}. + Further arguments that will be passed to \code{\link{mkinfit}}. } } \value{ diff --git a/man/plot.mkinfit.Rd b/man/plot.mkinfit.Rd index 89e5fb3c..a54bbbf6 100644 --- a/man/plot.mkinfit.Rd +++ b/man/plot.mkinfit.Rd @@ -6,33 +6,33 @@ } \description{ Solves the differential equations with the optimised and fixed parameters - from a previous successful call to \code{\link{mkinfit}} and plots + from a previous successful call to \code{\link{mkinfit}} and plots the observed data together with the solution of the fitted model. } \usage{ \method{plot}{mkinfit}(x, fit = x, obs_vars = names(fit$mkinmod$map), - xlab = "Time", ylab = "Observed", + xlab = "Time", ylab = "Observed", xlim = range(fit$data$time), ylim = "default", - col_obs = 1:length(obs_vars), pch_obs = col_obs, + col_obs = 1:length(obs_vars), pch_obs = col_obs, lty_obs = rep(1, length(obs_vars)), add = FALSE, legend = !add, show_residuals = FALSE, maxabs = "auto", sep_obs = FALSE, rel.height.middle = 0.9, - lpos = "topright", inset = c(0.05, 0.05), + lpos = "topright", inset = c(0.05, 0.05), show_errmin = FALSE, errmin_digits = 3, \dots) plot_sep(fit, sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE, \dots) } \arguments{ \item{x}{ - Alias for fit introduced for compatibility with the generic S3 method. + Alias for fit introduced for compatibility with the generic S3 method. } \item{fit}{ An object of class \code{\link{mkinfit}}. } \item{obs_vars}{ - A character vector of names of the observed variables for which the + 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. } @@ -66,7 +66,7 @@ plot_sep(fit, sep_obs = TRUE, show_residuals = TRUE, show_errmin = TRUE, \dots) \item{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 + "sep_obs" is given, the residual plots will be located to the right of the plots of the fitted curves. } \item{maxabs}{ diff --git a/man/plot.mmkin.Rd b/man/plot.mmkin.Rd index cfc7a35b..57f69c90 100644 --- a/man/plot.mmkin.Rd +++ b/man/plot.mmkin.Rd @@ -5,11 +5,11 @@ } \description{ 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 + 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. } \usage{ -\method{plot}{mmkin}(x, main = "auto", legends = 1, errmin_var = "All data", errmin_digits = 3, +\method{plot}{mmkin}(x, main = "auto", legends = 1, errmin_var = "All data", errmin_digits = 3, cex = 0.7, rel.height.middle = 0.9, ...) } \arguments{ @@ -46,12 +46,12 @@ } \examples{ # Only use one core not to offend CRAN checks, use Levenberg-Marquardt for speed - fits <- mmkin(c("FOMC", "HS"), list("FOCUS B" = FOCUS_2006_B, "FOCUS C" = FOCUS_2006_C), + fits <- mmkin(c("FOMC", "HS"), list("FOCUS B" = FOCUS_2006_B, "FOCUS C" = FOCUS_2006_C), cores = 1, quiet = TRUE, method.modFit = "Marq") plot(fits[, "FOCUS C"]) plot(fits["FOMC", ]) - # We can also plot a single fit, if we like the way mmkin works, but then the plot + # We can also plot a single fit, if we like the way mmkin works, but then the plot # height should be smaller than the plot width (this is not possible for the html pages # generated by staticdocs, as far as I know). plot(fits["FOMC", "FOCUS C"]) # same as plot(fits[1, 2]) diff --git a/man/transform_odeparms.Rd b/man/transform_odeparms.Rd index c6034424..71c179d7 100644 --- a/man/transform_odeparms.Rd +++ b/man/transform_odeparms.Rd @@ -5,7 +5,7 @@ Functions to transform and backtransform kinetic parameters for fitting } \description{ - The transformations are intended to map parameters that should only take + 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 @@ -17,9 +17,9 @@ This is no problem for the internal use in \code{\link{mkinfit}}. } \usage{ -transform_odeparms(parms, mkinmod, +transform_odeparms(parms, mkinmod, transform_rates = TRUE, transform_fractions = TRUE) -backtransform_odeparms(transparms, mkinmod, +backtransform_odeparms(transparms, mkinmod, transform_rates = TRUE, transform_fractions = TRUE) } \arguments{ @@ -32,22 +32,22 @@ backtransform_odeparms(transparms, mkinmod, \item{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. + before \code{\link{ilr}} transformation, the parameter names and + the information if the pathway to sink is included in the model. } \item{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 + 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. } \item{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 + 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. @@ -70,7 +70,7 @@ summary(fit, data=FALSE) # See transformed and backtransformed parameters \dontrun{ fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE) -summary(fit.2, data=FALSE) +summary(fit.2, data=FALSE) } initials <- fit$start$value diff --git a/test.R b/test.R deleted file mode 100644 index 1a81652d..00000000 --- a/test.R +++ /dev/null @@ -1,45 +0,0 @@ -library(pkgdown) -library(testthat) - -debug(pkgdown:::topic_index) - -build_reference_index(".") - -extract_tag(rd[[1]], "tag_title") # gives "\n" -extract_tag(rd[[4]], "tag_title") # gives the title - -extract_tag(rd[[1]], "tag_alias") # one aliase -extract_tag(rd[[13]], "tag_alias") # two aliases - -library(magrittr) - -extract_tag_1 <- function(x, tag) { - x %>% - purrr::keep(inherits, tag) %>% - unlist %>% - paste(collapse = "") %>% - trimws -} - -extract_tag <- function(x, tag) { - x %>% - purrr::keep(inherits, tag) %>% - purrr::map_chr(c(1, 1)) -} - -extract_tag_1 <- function(x, tag) { - x %>% - purrr::keep(inherits, tag) %>% - purrr::map_chr(function(x) trimws(paste(x, collapse = " "))) -} - -subset(rd[[1]][[1]], - -sapply(extract_tag_1(rd[[1]], "tag_title"), function(x) trimws(paste(x, collapse = " "))) - -extract_tag_1(rd[[1]], "tag_title") -extract_tag_1(rd[[4]], "tag_title") - -extract_tag_1(rd[[1]], "tag_alias") # one aliase -extract_tag_1(rd[[13]], "tag_alias") # two aliases - diff --git a/vignettes/FOCUS_D.Rmd b/vignettes/FOCUS_D.Rmd index 998cefa4..40bcb157 100644 --- a/vignettes/FOCUS_D.Rmd +++ b/vignettes/FOCUS_D.Rmd @@ -18,7 +18,7 @@ library(knitr) opts_chunk$set(tidy = FALSE, cache = TRUE)
```
-This is just a very simple vignette showing how to fit a degradation model for a parent
+This is just a very simple vignette showing how to fit a degradation model for a parent
compound with one transformation product using `mkin`. After loading the
library we look a the data. We have observed concentrations in the column named
`value` at the times specified in column `time` for the two observed variables
@@ -33,7 +33,7 @@ print(FOCUS_2006_D) Next we specify the degradation model: The parent compound degrades with simple first-order
kinetics (SFO) to one metabolite named m1, which also degrades with SFO kinetics.
-The call to mkinmod returns a degradation model. The differential equations represented in
+The call to mkinmod returns a degradation model. The differential equations represented in
R code can be found in the character vector `$diffs` of the `mkinmod` object. If
a C compiler (gcc) is installed and functional, the differential equation model will
be compiled from auto-generated C code.
@@ -51,7 +51,7 @@ fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) ```
A plot of the fit including a residual plot for both observed variables is obtained
-using the `plot_sep` method for `mkinfit` objects, which shows separate graphs for
+using the `plot_sep` method for `mkinfit` objects, which shows separate graphs for
all compounds and their residuals.
```{r plot, fig.height = 6, fig.width = 8}
diff --git a/vignettes/FOCUS_Z.Rnw b/vignettes/FOCUS_Z.Rnw index 5a66db24..5abda0e1 100644 --- a/vignettes/FOCUS_Z.Rnw +++ b/vignettes/FOCUS_Z.Rnw @@ -3,7 +3,7 @@ \documentclass[12pt,a4paper]{article} \usepackage{a4wide} \input{header} -\hypersetup{ +\hypersetup{ pdftitle = {Example evaluation of FOCUS dataset Z}, pdfsubject = {Manuscript}, pdfauthor = {Johannes Ranke}, @@ -50,15 +50,15 @@ report \citep{FOCUSkinetics2011}, p.350. require(mkin) LOD = 0.5 FOCUS_2006_Z = data.frame( - t = c(0, 0.04, 0.125, 0.29, 0.54, 1, 2, 3, 4, 7, 10, 14, 21, + t = c(0, 0.04, 0.125, 0.29, 0.54, 1, 2, 3, 4, 7, 10, 14, 21, 42, 61, 96, 124), - Z0 = c(100, 81.7, 70.4, 51.1, 41.2, 6.6, 4.6, 3.9, 4.6, 4.3, 6.8, + Z0 = c(100, 81.7, 70.4, 51.1, 41.2, 6.6, 4.6, 3.9, 4.6, 4.3, 6.8, 2.9, 3.5, 5.3, 4.4, 1.2, 0.7), - Z1 = c(0, 18.3, 29.6, 46.3, 55.1, 65.7, 39.1, 36, 15.3, 5.6, 1.1, + Z1 = c(0, 18.3, 29.6, 46.3, 55.1, 65.7, 39.1, 36, 15.3, 5.6, 1.1, 1.6, 0.6, 0.5 * LOD, NA, NA, NA), - Z2 = c(0, NA, 0.5 * LOD, 2.6, 3.8, 15.3, 37.2, 31.7, 35.6, 14.5, + Z2 = c(0, NA, 0.5 * LOD, 2.6, 3.8, 15.3, 37.2, 31.7, 35.6, 14.5, 0.8, 2.1, 1.9, 0.5 * LOD, NA, NA, NA), - Z3 = c(0, NA, NA, NA, NA, 0.5 * LOD, 9.2, 13.1, 22.3, 28.4, 32.5, + Z3 = c(0, NA, NA, NA, NA, 0.5 * LOD, 9.2, 13.1, 22.3, 28.4, 32.5, 25.2, 17.2, 4.8, 4.5, 2.8, 4.4)) FOCUS_2006_Z_mkin <- mkin_wide_to_long(FOCUS_2006_Z) @@ -66,9 +66,9 @@ FOCUS_2006_Z_mkin <- mkin_wide_to_long(FOCUS_2006_Z) \section{Parent compound and one metabolite} -The next step is to set up the models used for the kinetic analysis. As the +The next step is to set up the models used for the kinetic analysis. As the simultaneous fit of parent and the first metabolite is usually straightforward, -Step 1 (SFO for parent only) is skipped here. We start with the model 2a, +Step 1 (SFO for parent only) is skipped here. We start with the model 2a, with formation and decline of metabolite Z1 and the pathway from parent directly to sink included (default in mkin). @@ -80,15 +80,15 @@ plot_sep(m.Z.2a) summary(m.Z.2a, data = FALSE)$bpar @ -As obvious from the parameter summary (the \texttt{bpar} component of the +As obvious from the parameter summary (the \texttt{bpar} component of the summary), the kinetic rate constant from parent compound Z to sink -is negligible. Accordingly, the exact magnitude of the fitted parameter +is negligible. Accordingly, the exact magnitude of the fitted parameter \texttt{log k\_Z0\_sink} is ill-defined and the covariance matrix is not -returned (not shown, would be visible in the complete summary). +returned (not shown, would be visible in the complete summary). This suggests, in agreement with the analysis in the FOCUS kinetics report, to simplify the model by removing the pathway to sink. -A similar result can be obtained when formation fractions are used in the model +A similar result can be obtained when formation fractions are used in the model formulation: <<FOCUS_2006_Z_fits_2, echo=TRUE, fig.height=6>>= @@ -104,16 +104,16 @@ summary(m.Z.2a.ff, data = FALSE)$bpar Here, the ilr transformed formation fraction fitted in the model takes a very large value, and the backtransformed formation fraction from parent Z to Z1 is practically unity. Again, the covariance matrix is not returned as the model is -overparameterised. +overparameterised. The simplified model is obtained by setting the list component \texttt{sink} to \texttt{FALSE}.\footnote{If the model formulation without formation fractions is used, the same effect can be obtained by fixing the parameter \texttt{k\_Z\_sink} -to a value of zero.} +to a value of zero.} In the following, we use the parameterisation with formation fractions in order -to be able to compare with the results in the FOCUS guidance, and as it -makes it easier to use parameters obtained in a previous fit when adding a further +to be able to compare with the results in the FOCUS guidance, and as it +makes it easier to use parameters obtained in a previous fit when adding a further metabolite. <<FOCUS_2006_Z_fits_3, echo=TRUE, fig.height=6>>= @@ -130,10 +130,10 @@ to sink, the formation fraction is internally fixed to unity. \section{Including metabolites Z2 and Z3} As suggested in the FOCUS report, the pathway to sink was removed for metabolite Z1 as -well in the next step. While this step appears questionable on the basis of the above results, it -is followed here for the purpose of comparison. Also, in the FOCUS report, it is +well in the next step. While this step appears questionable on the basis of the above results, it +is followed here for the purpose of comparison. Also, in the FOCUS report, it is assumed that there is additional empirical evidence that Z1 quickly and exclusively -hydrolyses to Z2. +hydrolyses to Z2. <<FOCUS_2006_Z_fits_5, echo=TRUE, fig.height=7>>= Z.5 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), @@ -151,9 +151,9 @@ accelerate the optimization. Z.FOCUS <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE), Z1 = mkinsub("SFO", "Z2", sink = FALSE), Z2 = mkinsub("SFO", "Z3"), - Z3 = mkinsub("SFO"), + Z3 = mkinsub("SFO"), use_of_ff = "max") -m.Z.FOCUS <- mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, +m.Z.FOCUS <- mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, parms.ini = m.Z.5$bparms.ode, quiet = TRUE) plot_sep(m.Z.FOCUS) @@ -167,13 +167,13 @@ transformed parameters, however. \section{Using the SFORB model for parent and metabolites} -As the FOCUS report states, there is a certain tailing of the time course of metabolite -Z3. Also, the time course of the parent compound is not fitted very well using the +As the FOCUS report states, there is a certain tailing of the time course of metabolite +Z3. Also, the time course of the parent compound is not fitted very well using the SFO model, as residues at a certain low level remain. -Therefore, an additional model is offered here, using the single first-order +Therefore, an additional model is offered here, using the single first-order reversible binding (SFORB) model for metabolite Z3. As expected, the $\chi^2$ -error level is lower for metabolite Z3 using this model and the graphical +error level is lower for metabolite Z3 using this model and the graphical fit for Z3 is improved. However, the covariance matrix is not returned. <<FOCUS_2006_Z_fits_7, echo=TRUE, fig.height=8>>= @@ -198,7 +198,7 @@ m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE) plot_sep(m.Z.mkin.3) @ -This results in a much better representation of the behaviour of the parent +This results in a much better representation of the behaviour of the parent compound Z0. Finally, Z3 is added as well. These models appear overparameterised (no @@ -209,7 +209,7 @@ Z.mkin.4 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), Z1 = mkinsub("SFO", "Z2", sink = FALSE), Z2 = mkinsub("SFO", "Z3"), Z3 = mkinsub("SFO")) -m.Z.mkin.4 <- mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, +m.Z.mkin.4 <- mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, parms.ini = m.Z.mkin.3$bparms.ode, quiet = TRUE) plot_sep(m.Z.mkin.4) @@ -224,7 +224,7 @@ Z.mkin.5 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE), Z1 = mkinsub("SFO", "Z2", sink = FALSE), Z2 = mkinsub("SFO", "Z3"), Z3 = mkinsub("SFORB")) -m.Z.mkin.5 <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, +m.Z.mkin.5 <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, parms.ini = m.Z.mkin.4$bparms.ode[1:4], quiet = TRUE) plot_sep(m.Z.mkin.5) @@ -233,10 +233,10 @@ plot_sep(m.Z.mkin.5) The summary view of the backtransformed parameters shows that we get no confidence intervals due to overparameterisation. As the optimized \texttt{k\_Z3\_bound\_free} is excessively small, it seems reasonable to fix it to -zero. +zero. <<FOCUS_2006_Z_fits_11a, echo=TRUE>>= -m.Z.mkin.5a <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, +m.Z.mkin.5a <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, parms.ini = c(m.Z.mkin.5$bparms.ode[1:7], k_Z3_bound_free = 0), fixed_parms = "k_Z3_bound_free", @@ -244,7 +244,7 @@ m.Z.mkin.5a <- mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, plot_sep(m.Z.mkin.5a) @ -As expected, the residual plots for Z0 and Z3 are more random than in the case of the +As expected, the residual plots for Z0 and Z3 are more random than in the case of the all SFO model for which they were shown above. In conclusion, the model \texttt{Z.mkin.5a} is proposed as the best-fit model for the dataset from Appendix 7 of the FOCUS report. @@ -261,8 +261,8 @@ The endpoints obtained with this model are endpoints(m.Z.mkin.5a) @ -It is clear the degradation rate of Z3 towards the end of the experiment -is very low as DT50\_Z3\_b2 (the second Eigenvalue of the system of two differential +It is clear the degradation rate of Z3 towards the end of the experiment +is very low as DT50\_Z3\_b2 (the second Eigenvalue of the system of two differential equations representing the SFORB system for Z3, corresponding to the slower rate constant of the DFOP model) is reported to be infinity. However, this appears to be a feature of the data. diff --git a/vignettes/compiled_models.Rmd b/vignettes/compiled_models.Rmd index 9fd39d81..18e1a462 100644 --- a/vignettes/compiled_models.Rmd +++ b/vignettes/compiled_models.Rmd @@ -25,7 +25,7 @@ This evaluation is taken from the example section of mkinfit. When using an mkin equal to or greater than 0.9-36 and a C compiler (gcc) is available, you will see
a message that the model is being compiled from autogenerated C code when
defining a model using mkinmod. The `mkinmod()` function checks for presence of
-the gcc compiler using
+the gcc compiler using
```{r check_gcc}
Sys.which("gcc")
@@ -48,12 +48,12 @@ the microbenchmark package. library("microbenchmark")
library("ggplot2")
mb.1 <- microbenchmark(
- "deSolve, not compiled" = mkinfit(SFO_SFO, FOCUS_2006_D,
- solution_type = "deSolve",
+ "deSolve, not compiled" = mkinfit(SFO_SFO, FOCUS_2006_D,
+ solution_type = "deSolve",
use_compiled = FALSE, quiet = TRUE),
- "Eigenvalue based" = mkinfit(SFO_SFO, FOCUS_2006_D,
+ "Eigenvalue based" = mkinfit(SFO_SFO, FOCUS_2006_D,
solution_type = "eigen", quiet = TRUE),
- "deSolve, compiled" = mkinfit(SFO_SFO, FOCUS_2006_D,
+ "deSolve, compiled" = mkinfit(SFO_SFO, FOCUS_2006_D,
solution_type = "deSolve", quiet = TRUE),
times = 3, control = list(warmup = 0))
@@ -62,8 +62,8 @@ print(mb.1) autoplot(mb.1)
```
-We see that using the compiled model is by a factor of
-`r round(smb.1[1, "median"]/smb.1[3, "median"], 1)`
+We see that using the compiled model is by a factor of
+`r round(smb.1[1, "median"]/smb.1[3, "median"], 1)`
faster than using the R version with the default ode solver, and it is even
faster than the Eigenvalue based solution implemented in R which does not need
iterative solution of the ODEs:
@@ -75,7 +75,7 @@ smb.1["median"]/smb.1["deSolve, compiled", "median"] ## Model that can not be solved with Eigenvalues
-This evaluation is also taken from the example section of mkinfit.
+This evaluation is also taken from the example section of mkinfit.
```{r benchmark_FOMC_SFO, fig.height = 3}
FOMC_SFO <- mkinmod(
@@ -83,7 +83,7 @@ FOMC_SFO <- mkinmod( m1 = mkinsub( "SFO"))
mb.2 <- microbenchmark(
- "deSolve, not compiled" = mkinfit(FOMC_SFO, FOCUS_2006_D,
+ "deSolve, not compiled" = mkinfit(FOMC_SFO, FOCUS_2006_D,
use_compiled = FALSE, quiet = TRUE),
"deSolve, compiled" = mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE),
times = 3, control = list(warmup = 0))
@@ -93,8 +93,8 @@ smb.2["median"]/smb.2["deSolve, compiled", "median"] autoplot(mb.2)
```
-Here we get a performance benefit of a factor of
-`r round(smb.2[1, "median"]/smb.2[2, "median"], 1)`
+Here we get a performance benefit of a factor of
+`r round(smb.2[1, "median"]/smb.2[2, "median"], 1)`
using the version of the differential equation model compiled from C code!
This vignette was built with mkin `r packageVersion("mkin")` on
diff --git a/vignettes/mkin.Rmd b/vignettes/mkin.Rmd index a3982df9..062bfdac 100644 --- a/vignettes/mkin.Rmd +++ b/vignettes/mkin.Rmd @@ -26,8 +26,8 @@ opts_chunk$set(engine='R', tidy=FALSE) # Abstract In the regulatory evaluation of chemical substances like plant protection products (pesticides), biocides and other chemicals, degradation data play an -important role. For the evaluation of pesticide degradation experiments, -detailed guidance has been developed, based on nonlinear optimisation. +important role. For the evaluation of pesticide degradation experiments, +detailed guidance has been developed, based on nonlinear optimisation. The `R` add-on package `mkin` [@pkg:mkin] implements fitting some of the models recommended in this guidance from within R and calculates some statistical measures for data series within one or more compartments, for parent and @@ -38,23 +38,23 @@ library(mkin) # Define the kinetic model m_SFO_SFO_SFO <- mkinmod(parent = mkinsub("SFO", "M1"), M1 = mkinsub("SFO", "M2"), - M2 = mkinsub("SFO"), + M2 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) # Produce model predictions using some arbitrary parameters sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) d_SFO_SFO_SFO <- mkinpredict(m_SFO_SFO_SFO, - c(k_parent = 0.03, - f_parent_to_M1 = 0.5, k_M1 = log(2)/100, + c(k_parent = 0.03, + f_parent_to_M1 = 0.5, k_M1 = log(2)/100, f_M1_to_M2 = 0.9, k_M2 = log(2)/50), c(parent = 100, M1 = 0, M2 = 0), sampling_times) -# Generate a dataset by adding normally distributed errors with +# Generate a dataset by adding normally distributed errors with # standard deviation 3, for two replicates at each sampling time d_SFO_SFO_SFO_err <- add_err(d_SFO_SFO_SFO, reps = 2, - sdfunc = function(x) 3, + sdfunc = function(x) 3, n = 1, seed = 123456789 ) # Fit the model to the dataset @@ -67,7 +67,7 @@ plot_sep(f_SFO_SFO_SFO, lpos = c("topright", "bottomright", "bottomright")) # Background Many approaches are possible regarding the evaluation of chemical degradation -data. +data. The now deprecated `kinfit` package [@pkg:kinfit] in `R` [@rcore2016] implements the approach recommended in the kinetics report provided by the @@ -91,10 +91,10 @@ models based on differential equations to data. The code was first uploaded to the BerliOS platform. When this was taken down, the version control history was imported into the R-Forge site, where the code -is still mirrored today (see *e.g.* +is still mirrored today (see *e.g.* [the initial commit on 11 May 2010](http://cgit.jrwb.de/mkin/commit/?id=30cbb4092f6d2d3beff5800603374a0d009ad770)). -At that time, the R package `FME` (Flexible Modelling Environment) +At that time, the R package `FME` (Flexible Modelling Environment) [@soetaert2010] was already available, and provided a good basis for developing a package specifically tailored to the task. The remaining challenge was to make it as easy as possible for the users (including the author of this @@ -133,16 +133,16 @@ but until 2014, only simple first-order models could be specified for transformation products. Starting with KinGUII version 2.1, biphasic modelling of metabolites was also available in KinGUII. -A further graphical user interface (GUI) that has recently been brought to a decent -degree of maturity is the browser based GUI named `gmkin`. Please see its -[documentation page](http://kinfit.r-forge.r-project.org/gmkin_static) and +A further graphical user interface (GUI) that has recently been brought to a decent +degree of maturity is the browser based GUI named `gmkin`. Please see its +[documentation page](http://kinfit.r-forge.r-project.org/gmkin_static) and [manual](http://kinfit.r-forge.r-project.org/gmkin_static/vignettes/gmkin_manual.html) for further information. ## Recent developments Currently (June 2016), the main features available in `mkin` which are -not present in KinGUII or CAKE, are the speed increase by using compiled code when +not present in KinGUII or CAKE, are the speed increase by using compiled code when a compiler is present, parallel model fitting on multicore machines using the `mmkin` function, and the estimation of parameter confidence intervals based on transformed parameters. These are explained in more detail below. @@ -156,9 +156,9 @@ constants [compare @bates1988, p. 135], *i.e.* for their logarithms. Confidence intervals for the rate constants are then obtained using the appropriate backtransformation using the exponential function. -In the first version of `mkin` allowing for specifying models using +In the first version of `mkin` allowing for specifying models using formation fractions, a home-made reparameterisation was used in order to ensure -that the sum of formation fractions would not exceed unity. +that the sum of formation fractions would not exceed unity. This method is still used in the current version of KinGUII (v2.1 from April 2014), with a modification that allows for fixing the pathway to sink to zero. @@ -175,7 +175,7 @@ confidence intervals. ## Confidence intervals based on transformed parameters In the first attempt at providing improved parameter confidence intervals -introduced to `mkin` in 2013, confidence intervals obtained from +introduced to `mkin` in 2013, confidence intervals obtained from FME on the transformed parameters were simply all backtransformed one by one to yield asymetric confidence intervals for the backtransformed parameters. @@ -186,14 +186,14 @@ fractions that quantify the paths to each of the compounds formed from a specific parent compound, and no such 1:1 relation exists. Therefore, parameter confidence intervals for formation fractions obtained with -this method only appear valid for the case of a single transformation product, where +this method only appear valid for the case of a single transformation product, where only one formation fraction is to be estimated, directly corresponding to one component of the ilr transformed parameter. -The confidence intervals obtained by backtransformation for the cases where a -1:1 relation between transformed and original parameter exist are considered by +The confidence intervals obtained by backtransformation for the cases where a +1:1 relation between transformed and original parameter exist are considered by the author of this vignette to be more accurate than those obtained using a -re-estimation of the Hessian matrix after backtransformation, as implemented +re-estimation of the Hessian matrix after backtransformation, as implemented in the FME package. ## Parameter t-test based on untransformed parameters @@ -208,7 +208,7 @@ of the estimator for the parameters, is not fulfilled in the case of nonlinear r [@ranke2015]. However, this test is commonly used by industry, consultants and national authorities in order to decide on the reliability of parameter estimates, based on the FOCUS guidance mentioned above. Therefore, the results of this one-sided -t-test are included in the summary output from `mkin`. +t-test are included in the summary output from `mkin`. As it is not reasonable to test for significant difference of the transformed parameters (*e.g.* $log(k)$) from zero, the t-test is calculated based on the |