aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/SFORB.solution.R4
-rw-r--r--R/add_err.R6
-rw-r--r--R/endpoints.R12
-rw-r--r--R/ilr.R4
-rw-r--r--R/mkin_long_to_wide.R4
-rw-r--r--R/mkinds.R2
-rw-r--r--R/mkinerrmin.R14
-rw-r--r--R/mkinfit.R90
-rw-r--r--R/mkinmod.R20
-rw-r--r--R/mkinparplot.R22
-rw-r--r--R/mkinpredict.R26
-rw-r--r--R/mkinresplot.R10
-rw-r--r--R/mmkin.R10
-rw-r--r--R/plot.mkinfit.R32
-rw-r--r--R/plot.mmkin.R10
-rw-r--r--R/transform_odeparms.R10
16 files changed, 136 insertions, 140 deletions
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
}
diff --git a/R/ilr.R b/R/ilr.R
index 389653e7..620afc49 100644
--- a/R/ilr.R
+++ b/R/ilr.R
@@ -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
diff --git a/R/mkinds.R b/R/mkinds.R
index cb3c712a..60924c21 100644
--- a/R/mkinds.R
+++ b/R/mkinds.R
@@ -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])
}
}
diff --git a/R/mmkin.R b/R/mmkin.R
index be4526af..63542f4f 100644
--- a/R/mmkin.R
+++ b/R/mmkin.R
@@ -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

Contact - Imprint