From d2c1ab854491ff047135fa8377400a68499e72de Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 17 Jul 2014 12:53:30 +0200 Subject: Handle non-convergence and maximum number of iterations For details see NEWS.md --- R/mkinfit.R | 117 ++++++++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 82 insertions(+), 35 deletions(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index c6e13b97..d591c42a 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -28,7 +28,8 @@ mkinfit <- function(mkinmod, observed, fixed_initials = names(mkinmod$diffs)[-1], solution_type = "auto", method.ode = "lsoda", - method.modFit = "Marq", + method.modFit = c("Marq", "Port", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"), + maxit.modFit = "auto", control.modFit = list(), transform_rates = TRUE, transform_fractions = TRUE, @@ -40,6 +41,16 @@ mkinfit <- function(mkinmod, observed, trace_parms = FALSE, ...) { + # Check optimisation method and set maximum number of iterations if specified + method.modFit = match.arg(method.modFit) + if (maxit.modFit != "auto") { + if (method.modFit == "Marq") control.modFit$maxiter = maxit.modFit + if (method.modFit == "Port") control.modFit$iter.max = maxit.modFit + if (method.modFit %in% c("SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) { + control.modFit$maxit = maxit.modFit + } + } + # Get the names of the state variables in the model mod_vars <- names(mkinmod$diffs) @@ -281,48 +292,76 @@ mkinfit <- function(mkinmod, observed, upper[other_fraction_parms] <- 1 } - 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) - # if requested by the user - weight.ini = weight - if (!is.null(err)) weight.ini = "manual" - - if (!is.null(reweight.method)) { - if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented") - if(!quiet) { - cat("IRLS based on variance estimates for each observed variable\n") - } - if (!quiet) { - cat("Initial variance estimates are:\n") - print(signif(fit$var_ms_unweighted, 8)) - } - reweight.diff = 1 - n.iter <- 0 - if (!is.null(err)) observed$err.ini <- observed[[err]] - err = "err.irls" - while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) { - n.iter <- n.iter + 1 - sigma.old <- sqrt(fit$var_ms_unweighted) - observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)] - fit <- modFit(cost, fit$par, method = method.modFit, - control = control.modFit, lower = lower, upper = upper, ...) - reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2) + # 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, + lower = lower, upper = upper, ...) + + # Reiterate the fit until convergence of the variance components (IRLS) + # if requested by the user + weight.ini = weight + if (!is.null(err)) weight.ini = "manual" + + if (!is.null(reweight.method)) { + if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented") + if(!quiet) { + cat("IRLS based on variance estimates for each observed variable\n") + } if (!quiet) { - cat("Iteration", n.iter, "yields variance estimates:\n") + cat("Initial variance estimates are:\n") print(signif(fit$var_ms_unweighted, 8)) - cat("Sum of squared differences to last variance estimates:", - signif(reweight.diff, 2), "\n") + } + reweight.diff = 1 + n.iter <- 0 + if (!is.null(err)) observed$err.ini <- observed[[err]] + err = "err.irls" + while (reweight.diff > reweight.tol & n.iter < reweight.max.iter) { + n.iter <- n.iter + 1 + sigma.old <- sqrt(fit$var_ms_unweighted) + observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)] + fit <- modFit(cost, fit$par, method = method.modFit, + control = control.modFit, lower = lower, upper = upper, ...) + reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2) + if (!quiet) { + cat("Iteration", n.iter, "yields variance estimates:\n") + print(signif(fit$var_ms_unweighted, 8)) + cat("Sum of squared differences to last variance estimates:", + signif(reweight.diff, 2), "\n") + } } } + }) + + # Check for convergence + if (method.modFit == "Marq") { + if (!fit$info %in% c(1, 2, 3)) { + fit$warning = paste0("Optimisation by method ", method.modFit, + " did not converge.\n", + "The message returned by nls.lm is:\n", + fit$message) + warning(fit$warning) + } + } + 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, + " did not converge.\n", + "Convergence code is ", fit$convergence, + ifelse(is.null(fit$message), "", + paste0("\nMessage is ", fit$message))) + warning(fit$warning) + } } # We need to return some more data for summary and plotting fit$solution_type <- solution_type fit$transform_rates <- transform_rates fit$transform_fractions <- transform_fractions + fit$method.modFit <- method.modFit + fit$maxit.modFit <- maxit.modFit + fit$calls <- calls + fit$time <- fit_time # We also need the model for summary and plotting fit$mkinmod <- mkinmod @@ -449,6 +488,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, 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, @@ -461,6 +502,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, cov.scaled = covar * resvar, info = object$info, niter = object$iterations, + calls = object$calls, + time = object$time, stopmess = message, par = param, bpar = bparam) @@ -491,13 +534,17 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("Date of fit: ", x$date.fit, "\n") cat("Date of summary:", x$date.summary, "\n") + if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n") + cat("\nEquations:\n") print(noquote(as.character(x[["diffs"]]))) df <- x$df rdf <- df[2] - cat("\nMethod used for solution of differential equation system:\n") - cat(x$solution_type, "\n") + cat("\nModel predictions using solution type", x$solution_type, "\n") + + cat("\nFitted with method", x$method.modFit, + "using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n") cat("\nWeighting:", x$weight.ini) if(!is.null(x$reweight.method)) cat(" then iterative reweighting method", -- cgit v1.2.1 From 8def5006fc81c032c3fc99751e062cdb32a81cc1 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 17 Jul 2014 14:34:20 +0200 Subject: Return complete list of initial states after fitting This is useful for specifying state.ini in a subsequent call to mkinfit --- R/mkinfit.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index d591c42a..a8fbfc78 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -417,8 +417,11 @@ mkinfit <- function(mkinmod, observed, fit$bparms.optim <- bparms.optim fit$bparms.fixed <- bparms.fixed - # Return ode parameters for further fitting + # Return ode and state parameters for further fitting fit$bparms.ode <- bparms.all[mkinmod$parms] + fit$bparms.state <- c(bparms.all[setdiff(names(bparms.all), names(fit$bparms.ode))], + state.ini.fixed) + names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state)) fit$date <- date() -- cgit v1.2.1 From a1567638a3ba9f4d62fa199525097a94ddfd7912 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 21 Jul 2014 08:20:44 +0200 Subject: Bugfix, model shorthand, state.ini[[1]] from observed data - The bug occurred when using transform_rates=FALSE for FOMC, DFOP or HS - Make it possible to use mkinfit("SFO", ...) - Take initial mean value at time zero for the variable with the highest value in the observed data - Update of vignette/FOCUS_L - Improve the Makefile to build single vignettes --- R/mkinfit.R | 27 ++++++++++++++++++++++++--- R/transform_odeparms.R | 22 +++++++++++++++------- 2 files changed, 39 insertions(+), 10 deletions(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index a8fbfc78..39d084cb 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -23,7 +23,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "value")) mkinfit <- function(mkinmod, observed, parms.ini = "auto", - state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), + state.ini = "auto", fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], solution_type = "auto", @@ -41,6 +41,21 @@ mkinfit <- function(mkinmod, observed, trace_parms = FALSE, ...) { + # Check mkinmod and generate a model for the variable whithe the highest value + # if a suitable string is given + parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB") + presumed_parent_name = observed[which.max(observed$value), "name"] + if (class(mkinmod) != "mkinmod") { + if (mkinmod[[1]] %in% parent_models_available) { + speclist <- list(list(type = mkinmod, sink = TRUE)) + names(speclist) <- presumed_parent_name + mkinmod <- mkinmod(speclist = speclist) + } else { + stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ", + paste(parent_models_available, collapse = ", ")) + } + } + # Check optimisation method and set maximum number of iterations if specified method.modFit = match.arg(method.modFit) if (maxit.modFit != "auto") { @@ -55,7 +70,7 @@ mkinfit <- function(mkinmod, observed, mod_vars <- names(mkinmod$diffs) # Get the names of observed variables - obs_vars = names(mkinmod$spec) + obs_vars <- names(mkinmod$spec) # Subset observed data with names of observed data in the model observed <- subset(observed, name %in% obs_vars) @@ -137,6 +152,12 @@ mkinfit <- function(mkinmod, observed, } } + # Set default for state.ini if appropriate + if (state.ini[1] == "auto") { + state.ini = c(mean(subset(observed, time == 0 & name == presumed_parent_name)$value), + rep(0, length(mkinmod$diffs) - 1)) + } + # Name the inital state variable values if they are not named yet if(is.null(names(state.ini))) names(state.ini) <- mod_vars @@ -279,7 +300,7 @@ mkinfit <- function(mkinmod, observed, if (!transform_rates) { index_k <- grep("^k_", names(lower)) lower[index_k] <- 0 - other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2"), names(lower)) + other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb"), names(lower)) lower[other_rate_parms] <- 0 } diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 912a5c0a..f518ae32 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -69,7 +69,11 @@ transform_odeparms <- function(parms, mkinmod, # and HS parameter tb if transformation of rates is requested for (pname in c("alpha", "beta", "k1", "k2", "tb")) { if (!is.na(parms[pname])) { - transparms[paste0("log_", pname)] <- ifelse(transform_rates, log(parms[pname]), parms[pname]) + if (transform_rates) { + transparms[paste0("log_", pname)] <- log(parms[pname]) + } else { + transparms[pname] <- parms[pname] + } } } if (!is.na(parms["g"])) { @@ -130,12 +134,16 @@ backtransform_odeparms <- function(transparms, mkinmod, # Transform parameters also for FOMC, DFOP and HS models for (pname in c("alpha", "beta", "k1", "k2", "tb")) { - pname_trans = paste0("log_", pname) - if (!is.na(transparms[pname_trans])) { - parms[pname] <- ifelse(transform_rates, - exp(transparms[pname_trans]), - transparms[pname]) - } + if (transform_rates) { + pname_trans = paste0("log_", pname) + if (!is.na(transparms[pname_trans])) { + parms[pname] <- exp(transparms[pname_trans]) + } + } else { + if (!is.na(transparms[pname])) { + parms[pname] <- transparms[pname] + } + } } if (!is.na(transparms["g_ilr"])) { g_ilr <- transparms["g_ilr"] -- cgit v1.2.1 From 4c6f29fe2a3ece5a85160b891c89ce0f55299c11 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 23 Jul 2014 08:34:59 +0200 Subject: Parallel metabolite formation with formation fractions in mkinerrmin --- R/mkinerrmin.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) (limited to 'R') diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 9ebac6a4..09724730 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -65,11 +65,12 @@ mkinerrmin <- function(fit, alpha = 0.05) # Formation fractions are attributed to the target variable, so look # for source compartments with formation fractions for (source_var in fit$obs_vars) { + n.ff.source = length(grep(paste("^f", source_var, sep = "_"), + names(parms.optim))) + n.paths.source = length(fit$mkinmod$spec[[source_var]]$to) for (target_var in fit$mkinmod$spec[[source_var]]$to) { if (obs_var == target_var) { - n.ff.optim <- n.ff.optim + - length(grep(paste("^f", source_var, sep = "_"), - names(parms.optim))) + n.ff.optim <- n.ff.optim + n.ff.source/n.paths.source } } } -- cgit v1.2.1 From 8be01094ecbc65d4bdf1e7f819ef01b7a252b39d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 24 Jul 2014 09:58:55 +0200 Subject: Avoid artificial zero residual in mkinresplot --- R/mkinresplot.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'R') diff --git a/R/mkinresplot.R b/R/mkinresplot.R index c9a801fd..82ffd2cb 100644 --- a/R/mkinresplot.R +++ b/R/mkinresplot.R @@ -36,7 +36,8 @@ mkinresplot <- function (object, col_obs <- pch_obs <- 1:length(obs_vars) names(col_obs) <- names(pch_obs) <- obs_vars - plot(0, xlab = xlab, ylab = ylab, + plot(0, type = "n", + xlab = xlab, ylab = ylab, xlim = xlim, ylim = c(-1.2 * maxabs, 1.2 * maxabs), ...) -- cgit v1.2.1 From 9b947f0358d3a1b1fc922bfd0187ca444ce5811d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 24 Jul 2014 14:42:59 +0200 Subject: Bump version, better default for state.ini --- R/mkinfit.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index 39d084cb..c98c7586 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -154,8 +154,14 @@ mkinfit <- function(mkinmod, observed, # Set default for state.ini if appropriate if (state.ini[1] == "auto") { - state.ini = c(mean(subset(observed, time == 0 & name == presumed_parent_name)$value), - rep(0, length(mkinmod$diffs) - 1)) + presumed_parent_time_0 = subset(observed, + time == 0 & name == presumed_parent_name)$value + presumed_parent_time_0_mean = mean(presumed_parent_time_0, na.rm = TRUE) + if (is.na(presumed_parent_time_0_mean)) { + state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)) + } else { + state.ini = c(presumed_parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1)) + } } # Name the inital state variable values if they are not named yet -- cgit v1.2.1 From e5c955f82adf6139d76f842a0b85e5d383685793 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 25 Jul 2014 09:36:15 +0200 Subject: Fix internal naming of g for transform_fractions=FALSE --- R/transform_odeparms.R | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) (limited to 'R') diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index f518ae32..e64bac59 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -76,9 +76,15 @@ transform_odeparms <- function(parms, mkinmod, } } } + + # DFOP parameter g is treated as a fraction if (!is.na(parms["g"])) { g <- parms["g"] - transparms["g_ilr"] <- ifelse(transform_fractions, ilr(c(g, 1 - g)), g) + if (transform_fractions) { + transparms["g_ilr"] <- ilr(c(g, 1 - g)) + } else { + transparms["g"] <- g + } } return(transparms) @@ -145,9 +151,14 @@ backtransform_odeparms <- function(transparms, mkinmod, } } } + + # DFOP parameter g is treated as a fraction if (!is.na(transparms["g_ilr"])) { g_ilr <- transparms["g_ilr"] - parms["g"] <- ifelse(transform_fractions, invilr(g_ilr)[1], g_ilr) + parms["g"] <- invilr(g_ilr)[1] + } + if (!is.na(transparms["g"])) { + parms["g"] <- transparms["g"] } return(parms) -- cgit v1.2.1 From 37a252cb44fed78c4f7a00a2f7874f1c47456468 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 19 Aug 2014 17:55:06 +0200 Subject: Improve formatting of differential equations in output Rebuild of FOCUS_Z vignette with improved formatting --- R/mkinfit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index c98c7586..92c8f8bf 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -567,7 +567,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n") cat("\nEquations:\n") - print(noquote(as.character(x[["diffs"]]))) + cat(noquote(strwrap(x[["diffs"]], exdent = 11)), fill = TRUE) df <- x$df rdf <- df[2] -- cgit v1.2.1 From 81dc12a8c149e11d57190ebf59f3f7e9c7b99af7 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 20 Aug 2014 14:42:20 +0200 Subject: Fix mkinfit for a special situation See NEWS.md: When the parent was not the (only) variable with the highest value out of all variables in the observed data, it could happen that the initial value (state.ini) was not initialised. --- R/mkinfit.R | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index 92c8f8bf..0f488ec8 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -41,11 +41,11 @@ mkinfit <- function(mkinmod, observed, trace_parms = FALSE, ...) { - # Check mkinmod and generate a model for the variable whithe the highest value + # 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") - presumed_parent_name = observed[which.max(observed$value), "name"] if (class(mkinmod) != "mkinmod") { + presumed_parent_name = observed[which.max(observed$value), "name"] if (mkinmod[[1]] %in% parent_models_available) { speclist <- list(list(type = mkinmod, sink = TRUE)) names(speclist) <- presumed_parent_name @@ -153,14 +153,14 @@ mkinfit <- function(mkinmod, observed, } # Set default for state.ini if appropriate + parent_name = names(mkinmod$spec)[[1]] if (state.ini[1] == "auto") { - presumed_parent_time_0 = subset(observed, - time == 0 & name == presumed_parent_name)$value - presumed_parent_time_0_mean = mean(presumed_parent_time_0, na.rm = TRUE) - if (is.na(presumed_parent_time_0_mean)) { + parent_time_0 = subset(observed, time == 0 & name == parent_name)$value + parent_time_0_mean = mean(parent_time_0, na.rm = TRUE) + if (is.na(parent_time_0_mean)) { state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)) } else { - state.ini = c(presumed_parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1)) + state.ini = c(parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1)) } } -- cgit v1.2.1 From f30472ecd2afea6bd2153b8ad2bb2f663f3a2742 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 25 Aug 2014 10:39:40 +0200 Subject: Bug fix and unit tests for mkinerrmin See NEWS.md for details --- R/mkinerrmin.R | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) (limited to 'R') diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 09724730..2697d0a0 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -36,10 +36,11 @@ mkinerrmin <- function(fit, alpha = 0.05) suffixes = c("_mean", "_pred")) errdata <- errdata[order(errdata$time, errdata$name), ] - # Any value that is set to exactly zero is not really an observed value - # Remove those at time 0 - those are caused by the FOCUS recommendation - # to set metabolites occurring at time 0 to 0 - errdata <- subset(errdata, !(time == 0 & value_mean == 0)) + # 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 + # 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)) n.optim.overall <- length(parms.optim) -- cgit v1.2.1 From 2dc4f2f87048c28252992fe70055e0d7652a61ef Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 29 Aug 2014 12:33:26 +0200 Subject: Avoid an unnecessary warning --- R/plot.mkinfit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R') diff --git a/R/plot.mkinfit.R b/R/plot.mkinfit.R index 4132d96c..31746fb8 100644 --- a/R/plot.mkinfit.R +++ b/R/plot.mkinfit.R @@ -31,7 +31,7 @@ plot.mkinfit <- function(x, fit = x, { if (add && show_residuals) stop("If adding to an existing plot we can not show residuals") - if (ylim == "default") { + if (ylim[[1]] == "default") { ylim = c(0, max(subset(fit$data, variable %in% obs_vars)$observed, na.rm = TRUE)) } -- cgit v1.2.1 From 45942cc099c258fa07f42691cf2e2ee163f3b6d6 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Mon, 6 Oct 2014 20:17:17 +0200 Subject: Avoid a warning produced by initial fits in gmkin --- R/mkinfit.R | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index 0f488ec8..5c1de846 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -594,13 +594,15 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("\nOptimised, transformed parameters:\n") print(signif(x$par, digits = digits)) - cat("\nParameter correlation:\n") - if (!is.null(x$cov.unscaled)){ - Corr <- cov2cor(x$cov.unscaled) - rownames(Corr) <- colnames(Corr) <- rownames(x$par) - print(Corr, digits = digits, ...) - } else { - cat("Could not estimate covariance matrix; singular system:\n") + if (x$niter != 0) { + cat("\nParameter correlation:\n") + if (!is.null(x$cov.unscaled)){ + Corr <- cov2cor(x$cov.unscaled) + rownames(Corr) <- colnames(Corr) <- rownames(x$par) + print(Corr, digits = digits, ...) + } else { + cat("Could not estimate covariance matrix; singular system:\n") + } } cat("\nResidual standard error:", -- cgit v1.2.1 From 4510a609159216041f10a33146534f5a8366ac76 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 14 Oct 2014 22:04:54 +0200 Subject: Further formatting improvement for differential equations --- R/mkinfit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index 5c1de846..a966cea6 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -567,7 +567,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n") cat("\nEquations:\n") - cat(noquote(strwrap(x[["diffs"]], exdent = 11)), fill = TRUE) + writeLines(strwrap(x[["diffs"]], exdent = 11)) df <- x$df rdf <- df[2] -- cgit v1.2.1 From 65d31e345f9e61e9d05584b24df6a01c6c6ed18d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 15 Oct 2014 01:13:48 +0200 Subject: Switch to using the Port algorithm per default --- R/mkinfit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index a966cea6..6494ea1e 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -28,7 +28,7 @@ mkinfit <- function(mkinmod, observed, fixed_initials = names(mkinmod$diffs)[-1], solution_type = "auto", method.ode = "lsoda", - method.modFit = c("Marq", "Port", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"), + method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"), maxit.modFit = "auto", control.modFit = list(), transform_rates = TRUE, -- cgit v1.2.1 From c0803a4e6660a6ad3a3a219476b6269bab528bfe Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 22 Oct 2014 09:18:41 +0200 Subject: Always include 0 on y axis when plotting during the fit --- R/mkinfit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index 6494ea1e..e7cab1ee 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -281,7 +281,7 @@ mkinfit <- function(mkinmod, observed, atol = atol, rtol = rtol, ...) plot(0, type="n", - xlim = range(observed$time), ylim = range(observed$value, na.rm=TRUE), + 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)) -- cgit v1.2.1 From 3516b626be1aeb639d0735e79449424d2e987d7a Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 12 Nov 2014 12:46:15 +0100 Subject: Fix a typo in the list of possible arguments Thanks to Michael Brauer for the hint --- R/mkinfit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R') diff --git a/R/mkinfit.R b/R/mkinfit.R index e7cab1ee..59598631 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -28,7 +28,7 @@ mkinfit <- function(mkinmod, observed, fixed_initials = names(mkinmod$diffs)[-1], solution_type = "auto", method.ode = "lsoda", - method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFSG", "CG", "L-BFGS-B"), + method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"), maxit.modFit = "auto", control.modFit = list(), transform_rates = TRUE, -- cgit v1.2.1