diff options
Diffstat (limited to 'R')
-rw-r--r-- | R/SFORB.solution.R | 4 | ||||
-rw-r--r-- | R/add_err.R | 6 | ||||
-rw-r--r-- | R/endpoints.R | 12 | ||||
-rw-r--r-- | R/ilr.R | 4 | ||||
-rw-r--r-- | R/mkin_long_to_wide.R | 4 | ||||
-rw-r--r-- | R/mkinds.R | 2 | ||||
-rw-r--r-- | R/mkinerrmin.R | 14 | ||||
-rw-r--r-- | R/mkinfit.R | 90 | ||||
-rw-r--r-- | R/mkinmod.R | 20 | ||||
-rw-r--r-- | R/mkinparplot.R | 22 | ||||
-rw-r--r-- | R/mkinpredict.R | 26 | ||||
-rw-r--r-- | R/mkinresplot.R | 10 | ||||
-rw-r--r-- | R/mmkin.R | 10 | ||||
-rw-r--r-- | R/plot.mkinfit.R | 32 | ||||
-rw-r--r-- | R/plot.mmkin.R | 10 | ||||
-rw-r--r-- | R/transform_odeparms.R | 10 |
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
}
@@ -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
|