aboutsummaryrefslogtreecommitdiff
path: root/R/mkinfit.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/mkinfit.R')
-rw-r--r--R/mkinfit.R1801
1 files changed, 897 insertions, 904 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R
index d182f5d0..17fd59d0 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -1,904 +1,897 @@
-# Copyright (C) 2010-2019 Johannes Ranke
-# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG
-# Contact: jranke@uni-bremen.de
-
-# This file is part of the R package mkin
-
-# mkin is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software
-# Foundation, either version 3 of the License, or (at your option) any later
-# version.
-
-# This program is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-# details.
-
-# You should have received a copy of the GNU General Public License along with
-# this program. If not, see <http://www.gnu.org/licenses/>
-if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
-
-mkinfit <- function(mkinmod, observed,
- parms.ini = "auto",
- state.ini = "auto",
- err.ini = "auto",
- fixed_parms = NULL,
- fixed_initials = names(mkinmod$diffs)[-1],
- from_max_mean = FALSE,
- solution_type = c("auto", "analytical", "eigen", "deSolve"),
- method.ode = "lsoda",
- use_compiled = "auto",
- control = list(eval.max = 300, iter.max = 200),
- transform_rates = TRUE,
- transform_fractions = TRUE,
- quiet = FALSE,
- atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
- error_model = c("const", "obs", "tc"),
- error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", "OLS"),
- reweight.tol = 1e-8, reweight.max.iter = 10,
- trace_parms = FALSE,
- ...)
-{
- # Check mkinmod and generate a model for the variable whith the highest value
- # if a suitable string is given
- parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")
- if (class(mkinmod) != "mkinmod") {
- presumed_parent_name = observed[which.max(observed$value), "name"]
- if (mkinmod[[1]] %in% parent_models_available) {
- speclist <- list(list(type = mkinmod, sink = TRUE))
- names(speclist) <- presumed_parent_name
- mkinmod <- mkinmod(speclist = speclist)
- } else {
- stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ",
- paste(parent_models_available, collapse = ", "))
- }
- }
-
- # Get the names of the state variables in the model
- mod_vars <- names(mkinmod$diffs)
-
- # Get the names of observed variables
- obs_vars <- names(mkinmod$spec)
-
- # Subset observed data with names of observed data in the model and remove NA values
- observed <- subset(observed, name %in% obs_vars)
- observed <- subset(observed, !is.na(value))
-
- # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model)
- if (any(observed$value == 0)) {
- warning("Observations with value of zero were removed from the data")
- observed <- subset(observed, value != 0)
- }
-
- # Obtain data for decline from maximum mean value if requested
- if (from_max_mean) {
- # This is only used for simple decline models
- if (length(obs_vars) > 1)
- stop("Decline from maximum is only implemented for models with a single observed variable")
-
- means <- aggregate(value ~ time, data = observed, mean, na.rm=TRUE)
- t_of_max <- means[which.max(means$value), "time"]
- observed <- subset(observed, time >= t_of_max)
- observed$time <- observed$time - t_of_max
- }
-
- # Number observations used for fitting
- n_observed <- nrow(observed)
-
- # Define starting values for parameters where not specified by the user
- if (parms.ini[[1]] == "auto") parms.ini = vector()
-
- # Warn for inital parameter specifications that are not in the model
- wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms)
- if (length(wrongpar.names) > 0) {
- warning("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "),
- " not used in the model")
- parms.ini <- parms.ini[setdiff(names(parms.ini), wrongpar.names)]
- }
-
- # Warn that the sum of formation fractions may exceed one if they are not
- # fitted in the transformed way
- if (mkinmod$use_of_ff == "max" & transform_fractions == FALSE) {
- warning("The sum of formation fractions may exceed one if you do not use ",
- "transform_fractions = TRUE." )
- for (box in mod_vars) {
- # Stop if formation fractions are not transformed and we have no sink
- if (mkinmod$spec[[box]]$sink == FALSE) {
- stop("If formation fractions are not transformed during the fitting, ",
- "it is not supported to turn off pathways to sink.\n ",
- "Consider turning on the transformation of formation fractions or ",
- "setting up a model with use_of_ff = 'min'.\n")
- }
- }
- }
-
- # Do not allow fixing formation fractions if we are using the ilr transformation,
- # this is not supported
- if (transform_fractions == TRUE && length(fixed_parms) > 0) {
- if (any(grepl("^f_", fixed_parms))) {
- stop("Fixing formation fractions is not supported when using the ilr ",
- "transformation.")
- }
- }
-
- # Set initial parameter values, including a small increment (salt)
- # to avoid linear dependencies (singular matrix) in Eigenvalue based solutions
- k_salt = 0
- defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
- for (parmname in defaultpar.names) {
- # Default values for rate constants, depending on the parameterisation
- if (grepl("^k", parmname)) {
- parms.ini[parmname] = 0.1 + k_salt
- k_salt = k_salt + 1e-4
- }
- # Default values for rate constants for reversible binding
- if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1
- if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02
- # Default values for IORE exponents
- if (grepl("^N", parmname)) parms.ini[parmname] = 1.1
- # Default values for the FOMC, DFOP and HS models
- if (parmname == "alpha") parms.ini[parmname] = 1
- if (parmname == "beta") parms.ini[parmname] = 10
- if (parmname == "k1") parms.ini[parmname] = 0.1
- if (parmname == "k2") parms.ini[parmname] = 0.01
- if (parmname == "tb") parms.ini[parmname] = 5
- if (parmname == "g") parms.ini[parmname] = 0.5
- if (parmname == "kmax") parms.ini[parmname] = 0.1
- if (parmname == "k0") parms.ini[parmname] = 0.0001
- if (parmname == "r") parms.ini[parmname] = 0.2
- }
- # Default values for formation fractions in case they are present
- for (box in mod_vars) {
- f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)]
- if (length(f_names) > 0) {
- # We need to differentiate between default and specified fractions
- # and set the unspecified to 1 - sum(specified)/n_unspecified
- f_default_names <- intersect(f_names, defaultpar.names)
- f_specified_names <- setdiff(f_names, defaultpar.names)
- sum_f_specified = sum(parms.ini[f_specified_names])
- if (sum_f_specified > 1) {
- stop("Starting values for the formation fractions originating from ",
- box, " sum up to more than 1.")
- }
- if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1
- else {
- n_unspecified = length(f_default_names)
- }
- parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified
- }
- }
-
- # Set default for state.ini if appropriate
- parent_name = names(mkinmod$spec)[[1]]
- if (state.ini[1] == "auto") {
- parent_time_0 = subset(observed, time == 0 & name == parent_name)$value
- parent_time_0_mean = mean(parent_time_0, na.rm = TRUE)
- if (is.na(parent_time_0_mean)) {
- state.ini = c(100, rep(0, length(mkinmod$diffs) - 1))
- } else {
- state.ini = c(parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1))
- }
- }
-
- # Name the inital state variable values if they are not named yet
- if(is.null(names(state.ini))) names(state.ini) <- mod_vars
-
- # Transform initial parameter values for fitting
- transparms.ini <- transform_odeparms(parms.ini, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Parameters to be optimised:
- # Kinetic parameters in parms.ini whose names are not in fixed_parms
- parms.fixed <- parms.ini[fixed_parms]
- parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)]
-
- transparms.fixed <- transform_odeparms(parms.fixed, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- transparms.optim <- transform_odeparms(parms.optim, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
-
- # Inital state variables in state.ini whose names are not in fixed_initials
- state.ini.fixed <- state.ini[fixed_initials]
- state.ini.optim <- state.ini[setdiff(names(state.ini), fixed_initials)]
-
- # Preserve names of state variables before renaming initial state variable
- # parameters
- state.ini.optim.boxnames <- names(state.ini.optim)
- state.ini.fixed.boxnames <- names(state.ini.fixed)
- if(length(state.ini.optim) > 0) {
- names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
- }
- if(length(state.ini.fixed) > 0) {
- names(state.ini.fixed) <- paste(names(state.ini.fixed), "0", sep="_")
- }
-
- # Decide if the solution of the model can be based on a simple analytical
- # formula, the spectral decomposition of the matrix (fundamental system)
- # or a numeric ode solver from the deSolve package
- # Prefer deSolve over eigen if a compiled model is present and use_compiled
- # is not set to FALSE
- solution_type = match.arg(solution_type)
- if (solution_type == "analytical" && length(mkinmod$spec) > 1)
- stop("Analytical solution not implemented for models with metabolites.")
- if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
- stop("Eigenvalue based solution not possible, coefficient matrix not present.")
- if (solution_type == "auto") {
- if (length(mkinmod$spec) == 1) {
- solution_type = "analytical"
- } else {
- if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) {
- solution_type = "deSolve"
- } else {
- if (is.matrix(mkinmod$coefmat)) {
- solution_type = "eigen"
- if (max(observed$value, na.rm = TRUE) < 0.1) {
- stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone")
- }
- } else {
- solution_type = "deSolve"
- }
- }
- }
- }
-
- # Get the error model and the algorithm for fitting
- err_mod <- match.arg(error_model)
- error_model_algorithm = match.arg(error_model_algorithm)
- if (error_model_algorithm == "OLS") {
- if (err_mod != "const") stop("OLS is only appropriate for constant variance")
- }
- if (error_model_algorithm == "auto") {
- error_model_algorithm = switch(err_mod,
- const = "OLS", obs = "d_3", tc = "d_3")
- }
- errparm_names <- switch(err_mod,
- "const" = "sigma",
- "obs" = paste0("sigma_", obs_vars),
- "tc" = c("sigma_low", "rsd_high"))
- errparm_names_optim <- if (error_model_algorithm == "OLS") NULL else errparm_names
-
- # Define starting values for the error model
- if (err.ini[1] != "auto") {
- if (!identical(names(err.ini), errparm_names)) {
- stop("Please supply initial values for error model components ", paste(errparm_names, collapse = ", "))
- } else {
- errparms = err.ini
- }
- } else {
- if (err_mod == "const") {
- errparms = 3
- }
- if (err_mod == "obs") {
- errparms = rep(3, length(obs_vars))
- }
- if (err_mod == "tc") {
- errparms <- c(sigma_low = 0.1, rsd_high = 0.1)
- }
- names(errparms) <- errparm_names
- }
- if (error_model_algorithm == "OLS") {
- errparms_optim <- NULL
- } else {
- errparms_optim <- errparms
- }
-
- # Define outtimes for model solution.
- # Include time points at which observed data are available
- outtimes = sort(unique(c(observed$time, seq(min(observed$time),
- max(observed$time),
- length.out = n.outtimes))))
-
- # Define the objective function for optimisation, including (back)transformations
- cost_function <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...)
- {
- assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter
-
- # Trace parameter values if requested and if we are actually optimising
- if(trace_parms & update_data) cat(P, "\n")
-
- # Determine local parameter values for the cost estimation
- if (is.numeric(fixed_degparms)) {
- cost_degparms <- fixed_degparms
- cost_errparms <- P
- degparms_fixed = TRUE
- } else {
- degparms_fixed = FALSE
- }
-
- if (is.numeric(fixed_errparms)) {
- cost_degparms <- P
- cost_errparms <- fixed_errparms
- errparms_fixed = TRUE
- } else {
- errparms_fixed = FALSE
- }
-
- if (OLS) {
- cost_degparms <- P
- cost_errparms <- numeric(0)
- }
-
- if (!OLS & !degparms_fixed & !errparms_fixed) {
- cost_degparms <- P[1:(length(P) - length(errparms))]
- cost_errparms <- P[(length(cost_degparms) + 1):length(P)]
- }
-
- # Initial states for t0
- if(length(state.ini.optim) > 0) {
- odeini <- c(cost_degparms[1:length(state.ini.optim)], state.ini.fixed)
- names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
- } else {
- odeini <- state.ini.fixed
- names(odeini) <- state.ini.fixed.boxnames
- }
-
- odeparms.optim <- cost_degparms[(length(state.ini.optim) + 1):length(cost_degparms)]
-
- if (trans == TRUE) {
- odeparms <- c(odeparms.optim, transparms.fixed)
- parms <- backtransform_odeparms(odeparms, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- } else {
- parms <- c(odeparms.optim, parms.fixed)
- }
-
- # Solve the system with current parameter values
- out <- mkinpredict(mkinmod, parms,
- odeini, outtimes,
- solution_type = solution_type,
- use_compiled = use_compiled,
- method.ode = method.ode,
- atol = atol, rtol = rtol, ...)
-
- out_long <- mkin_wide_to_long(out, time = "time")
-
- if (err_mod == "const") {
- observed$std <- if (OLS) NA else cost_errparms["sigma"]
- }
- if (err_mod == "obs") {
- std_names <- paste0("sigma_", observed$name)
- observed$std <- cost_errparms[std_names]
- }
- if (err_mod == "tc") {
- tmp <- merge(observed, out_long, by = c("time", "name"))
- tmp$name <- ordered(tmp$name, levels = obs_vars)
- tmp <- tmp[order(tmp$name, tmp$time), ]
- observed$std <- sqrt(cost_errparms["sigma_low"]^2 + tmp$value.y^2 * cost_errparms["rsd_high"]^2)
- }
-
- cost_data <- merge(observed[c("name", "time", "value", "std")], out_long,
- by = c("name", "time"), suffixes = c(".observed", ".predicted"))
-
- if (OLS) {
- # Cost is the sum of squared residuals
- cost <- with(cost_data, sum((value.observed - value.predicted)^2))
- } else {
- # Cost is the negative log-likelihood
- cost <- - with(cost_data,
- sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE)))
- }
-
- # We update the current cost and data during the optimisation, not
- # during hessian calculations
- if (update_data) {
-
- assign("out_predicted", out_long, inherits = TRUE)
- assign("current_data", cost_data, inherits = TRUE)
-
- if (cost < cost.current) {
- assign("cost.current", cost, inherits = TRUE)
- if (!quiet) cat(ifelse(OLS, "Sum of squared residuals", "Negative log-likelihood"),
- " at call ", calls, ": ", cost.current, "\n", sep = "")
- }
- }
- return(cost)
- }
-
- names_optim <- c(names(state.ini.optim),
- names(transparms.optim),
- errparm_names_optim)
- n_optim <- length(names_optim)
-
- # Define lower and upper bounds other than -Inf and Inf for parameters
- # for which no internal transformation is requested in the call to mkinfit
- # and for optimised error model parameters
- lower <- rep(-Inf, n_optim)
- upper <- rep(Inf, n_optim)
- names(lower) <- names(upper) <- names_optim
-
- # IORE exponents are not transformed, but need a lower bound
- index_N <- grep("^N", names(lower))
- lower[index_N] <- 0
-
- if (!transform_rates) {
- index_k <- grep("^k_", names(lower))
- lower[index_k] <- 0
- index_k__iore <- grep("^k__iore_", names(lower))
- lower[index_k__iore] <- 0
- other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb", "r"), names(lower))
- lower[other_rate_parms] <- 0
- }
-
- if (!transform_fractions) {
- index_f <- grep("^f_", names(upper))
- lower[index_f] <- 0
- upper[index_f] <- 1
- other_fraction_parms <- intersect(c("g"), names(upper))
- lower[other_fraction_parms] <- 0
- upper[other_fraction_parms] <- 1
- }
-
- if (err_mod == "const") {
- if (error_model_algorithm != "OLS") {
- lower["sigma"] <- 0
- }
- }
- if (err_mod == "obs") {
- index_sigma <- grep("^sigma_", names(lower))
- lower[index_sigma] <- 0
- }
- if (err_mod == "tc") {
- lower["sigma_low"] <- 0
- lower["rsd_high"] <- 0
- }
-
- # Counter for cost function evaluations
- calls = 0
- cost.current <- Inf
- out_predicted <- NA
- current_data <- NA
-
- # Show parameter names if tracing is requested
- if(trace_parms) cat(names_optim, "\n")
-
- # browser()
-
- # Do the fit and take the time until the hessians are calculated
- fit_time <- system.time({
- degparms <- c(state.ini.optim, transparms.optim)
- n_degparms <- length(degparms)
- degparms_index <- seq(1, n_degparms)
- errparms_index <- seq(n_degparms + 1, length.out = length(errparms))
-
- if (error_model_algorithm == "d_3") {
- if (!quiet) message("Directly optimising the complete model")
- parms.start <- c(degparms, errparms)
- fit_direct <- nlminb(parms.start, cost_function,
- lower = lower[names(parms.start)],
- upper = upper[names(parms.start)],
- control = control, ...)
- fit_direct$logLik <- - cost.current
- if (error_model_algorithm == "direct") {
- degparms <- fit_direct$par[degparms_index]
- errparms <- fit_direct$par[errparms_index]
- } else {
- cost.current <- Inf # reset to avoid conflict with the OLS step
- }
- }
- if (error_model_algorithm != "direct") {
- if (!quiet) message("Ordinary least squares optimisation")
- fit <- nlminb(degparms, cost_function, control = control,
- lower = lower[names(degparms)],
- upper = upper[names(degparms)], OLS = TRUE, ...)
- degparms <- fit$par
-
- # Get the maximum likelihood estimate for sigma at the optimum parameter values
- current_data$residual <- current_data$value.observed - current_data$value.predicted
- sigma_mle <- sqrt(sum(current_data$residual^2)/nrow(current_data))
-
- # Use that estimate for the constant variance, or as first guess if err_mod = "obs"
- if (err_mod != "tc") {
- errparms[names(errparms)] <- sigma_mle
- }
- fit$par <- c(fit$par, errparms)
-
- cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
- fit$logLik <- - cost.current
- }
- if (error_model_algorithm %in% c("threestep", "fourstep", "d_3")) {
- if (!quiet) message("Optimising the error model")
- fit <- nlminb(errparms, cost_function, control = control,
- lower = lower[names(errparms)],
- upper = upper[names(errparms)],
- fixed_degparms = degparms, ...)
- errparms <- fit$par
- }
- if (error_model_algorithm == "fourstep") {
- if (!quiet) message("Optimising the degradation model")
- fit <- nlminb(degparms, cost_function, control = control,
- lower = lower[names(degparms)],
- upper = upper[names(degparms)],
- fixed_errparms = errparms, ...)
- degparms <- fit$par
- }
- if (error_model_algorithm %in%
- c("direct", "twostep", "threestep", "fourstep", "d_3")) {
- if (!quiet) message("Optimising the complete model")
- parms.start <- c(degparms, errparms)
- fit <- nlminb(parms.start, cost_function,
- lower = lower[names(parms.start)],
- upper = upper[names(parms.start)],
- control = control, ...)
- degparms <- fit$par[degparms_index]
- errparms <- fit$par[errparms_index]
- fit$logLik <- - cost.current
-
- if (error_model_algorithm == "d_3") {
- d_3_messages = c(
- same = "Direct fitting and three-step fitting yield approximately the same likelihood",
- threestep = "Three-step fitting yielded a higher likelihood than direct fitting",
- direct = "Direct fitting yielded a higher likelihood than three-step fitting")
- rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
- if (rel_diff < 0.0001) {
- if (!quiet) message(d_3_messages["same"])
- fit$d_3_message <- d_3_messages["same"]
- } else {
- if (fit$logLik > fit_direct$logLik) {
- if (!quiet) message(d_3_messages["threestep"])
- fit$d_3_message <- d_3_messages["threestep"]
- } else {
- if (!quiet) message(d_3_messages["direct"])
- fit <- fit_direct
- fit$d_3_message <- d_3_messages["direct"]
- }
- }
- }
- }
- if (err_mod != "const" & error_model_algorithm == "IRLS") {
- reweight.diff <- 1
- n.iter <- 0
- errparms_last <- errparms
-
- while (reweight.diff > reweight.tol &
- n.iter < reweight.max.iter) {
-
- if (!quiet) message("Optimising the error model")
- fit <- nlminb(errparms, cost_function, control = control,
- lower = lower[names(errparms)],
- upper = upper[names(errparms)],
- fixed_degparms = degparms, ...)
- errparms <- fit$par
-
- if (!quiet) message("Optimising the degradation model")
- fit <- nlminb(degparms, cost_function, control = control,
- lower = lower[names(degparms)],
- upper = upper[names(degparms)],
- fixed_errparms = errparms, ...)
- degparms <- fit$par
-
- reweight.diff <- dist(rbind(errparms, errparms_last))
- errparms_last <- errparms
-
- fit$par <- c(fit$par, errparms)
- cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
- fit$logLik <- - cost.current
- }
- }
-
- fit$hessian <- try(numDeriv::hessian(cost_function, c(degparms, errparms), OLS = FALSE,
- update_data = FALSE), silent = TRUE)
-
- # Backtransform parameters
- bparms.optim = backtransform_odeparms(fit$par, mkinmod,
- transform_rates = transform_rates,
- transform_fractions = transform_fractions)
- bparms.fixed = c(state.ini.fixed, parms.fixed)
- bparms.all = c(bparms.optim, parms.fixed)
-
- fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.all, errparms),
- OLS = FALSE, trans = FALSE, update_data = FALSE), silent = TRUE)
- })
-
- fit$error_model_algorithm <- error_model_algorithm
-
- if (fit$convergence != 0) {
- fit$warning = paste0("Optimisation did not converge:\n", fit$message)
- warning(fit$warning)
- } else {
- if(!quiet) message("Optimisation successfully terminated.\n")
- }
-
- # We need to return some more data for summary and plotting
- fit$solution_type <- solution_type
- fit$transform_rates <- transform_rates
- fit$transform_fractions <- transform_fractions
- fit$reweight.tol <- reweight.tol
- fit$reweight.max.iter <- reweight.max.iter
- fit$control <- control
- fit$calls <- calls
- fit$time <- fit_time
-
- # We also need the model for summary and plotting
- fit$mkinmod <- mkinmod
-
- # We need data and predictions for summary and plotting
- fit$observed <- observed
- fit$obs_vars <- obs_vars
- fit$predicted <- out_predicted
-
- # Residual sum of squares as a function of the fitted parameters
- fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE)
-
- # Log-likelihood with possibility to fix degparms or errparms
- fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE) {
- - cost_function(P, fixed_degparms = fixed_degparms,
- fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE)
- }
-
- # Collect initial parameter values in three dataframes
- fit$start <- data.frame(value = c(state.ini.optim,
- parms.optim, errparms_optim))
- fit$start$type = c(rep("state", length(state.ini.optim)),
- rep("deparm", length(parms.optim)),
- rep("error", length(errparms_optim)))
-
- fit$start_transformed = data.frame(
- value = c(state.ini.optim, transparms.optim, errparms_optim),
- lower = lower,
- upper = upper)
-
- fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))
- fit$fixed$type = c(rep("state", length(state.ini.fixed)),
- rep("deparm", length(parms.fixed)))
-
- # Sort observed, predicted and residuals
- current_data$name <- ordered(current_data$name, levels = obs_vars)
-
- ordered_data <- current_data[order(current_data$name, current_data$time), ]
-
- fit$data <- data.frame(time = ordered_data$time,
- variable = ordered_data$name,
- observed = ordered_data$value.observed,
- predicted = ordered_data$value.predicted)
-
- fit$data$residual <- fit$data$observed - fit$data$predicted
-
- fit$atol <- atol
- fit$rtol <- rtol
- fit$err_mod <- err_mod
-
- # Return different sets of backtransformed parameters for summary and plotting
- fit$bparms.optim <- bparms.optim
- fit$bparms.fixed <- bparms.fixed
-
- # Return ode and state parameters for further fitting
- fit$bparms.ode <- bparms.all[mkinmod$parms]
- fit$bparms.state <- c(bparms.all[setdiff(names(bparms.all), names(fit$bparms.ode))],
- state.ini.fixed)
- names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state))
-
- fit$errparms <- errparms
- fit$df.residual <- n_observed - length(c(degparms, errparms))
-
- fit$date <- date()
- fit$version <- as.character(utils::packageVersion("mkin"))
- fit$Rversion <- paste(R.version$major, R.version$minor, sep=".")
-
- class(fit) <- c("mkinfit", "modFit")
- return(fit)
-}
-
-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
- param <- object$par
- pnames <- names(param)
- bpnames <- names(object$bparms.optim)
- epnames <- names(object$errparms)
- p <- length(param)
- mod_vars <- names(object$mkinmod$diffs)
- covar <- try(solve(object$hessian), silent = TRUE)
- covar_notrans <- try(solve(object$hessian_notrans), silent = TRUE)
- rdf <- object$df.residual
-
- if (!is.numeric(covar) | is.na(covar[1])) {
- covar <- NULL
- se <- lci <- uci <- rep(NA, p)
- } else {
- rownames(covar) <- colnames(covar) <- pnames
- se <- sqrt(diag(covar))
- lci <- param + qt(alpha/2, rdf) * se
- uci <- param + qt(1-alpha/2, rdf) * se
- }
-
- beparms.optim <- c(object$bparms.optim, object$par[epnames])
- if (!is.numeric(covar_notrans) | is.na(covar_notrans[1])) {
- covar_notrans <- NULL
- se_notrans <- tval <- pval <- rep(NA, p)
- } else {
- rownames(covar_notrans) <- colnames(covar_notrans) <- c(bpnames, epnames)
- se_notrans <- sqrt(diag(covar_notrans))
- tval <- beparms.optim / se_notrans
- pval <- pt(abs(tval), rdf, lower.tail = FALSE)
- }
-
- names(se) <- pnames
-
- param <- cbind(param, se, lci, uci)
- dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))
-
- bparam <- cbind(Estimate = beparms.optim, se_notrans,
- "t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA)
-
- # Transform boundaries of CI for one parameter at a time,
- # with the exception of sets of formation fractions (single fractions are OK).
- f_names_skip <- character(0)
- for (box in mod_vars) { # Figure out sets of fractions to skip
- f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
- n_paths <- length(f_names)
- if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
- }
-
- for (pname in pnames) {
- if (!pname %in% f_names_skip) {
- par.lower <- param[pname, "Lower"]
- par.upper <- param[pname, "Upper"]
- names(par.lower) <- names(par.upper) <- pname
- bpl <- backtransform_odeparms(par.lower, object$mkinmod,
- object$transform_rates,
- object$transform_fractions)
- bpu <- backtransform_odeparms(par.upper, object$mkinmod,
- object$transform_rates,
- object$transform_fractions)
- bparam[names(bpl), "Lower"] <- bpl
- bparam[names(bpu), "Upper"] <- bpu
- }
- }
- bparam[epnames, c("Lower", "Upper")] <- param[epnames, c("Lower", "Upper")]
-
- ans <- list(
- version = as.character(utils::packageVersion("mkin")),
- Rversion = paste(R.version$major, R.version$minor, sep="."),
- date.fit = object$date,
- date.summary = date(),
- solution_type = object$solution_type,
- warning = object$warning,
- use_of_ff = object$mkinmod$use_of_ff,
- error_model_algorithm = object$error_model_algorithm,
- df = c(p, rdf),
- covar = covar,
- covar_notrans = covar_notrans,
- err_mod = object$err_mod,
- niter = object$iterations,
- calls = object$calls,
- time = object$time,
- par = param,
- bpar = bparam)
-
- if (!is.null(object$version)) {
- ans$fit_version <- object$version
- ans$fit_Rversion <- object$Rversion
- }
-
- ans$diffs <- object$mkinmod$diffs
- if(data) ans$data <- object$data
- ans$start <- object$start
- ans$start_transformed <- object$start_transformed
-
- ans$fixed <- object$fixed
-
- ans$errmin <- mkinerrmin(object, alpha = 0.05)
-
- if (object$calls > 0) {
- if (!is.null(ans$covar)){
- Corr <- cov2cor(ans$covar)
- rownames(Corr) <- colnames(Corr) <- rownames(ans$par)
- ans$Corr <- Corr
- } else {
- warning("Could not calculate correlation; no covariance matrix")
- }
- }
-
- ans$bparms.ode <- object$bparms.ode
- ep <- endpoints(object)
- if (length(ep$ff) != 0)
- ans$ff <- ep$ff
- if (distimes) ans$distimes <- ep$distimes
- if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
- if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message
- class(ans) <- c("summary.mkinfit", "summary.modFit")
- return(ans)
-}
-
-# Expanded from print.summary.modFit
-print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), ...) {
- if (is.null(x$fit_version)) {
- cat("mkin version: ", x$version, "\n")
- cat("R version: ", x$Rversion, "\n")
- } else {
- cat("mkin version used for fitting: ", x$fit_version, "\n")
- cat("R version used for fitting: ", x$fit_Rversion, "\n")
- }
-
- cat("Date of fit: ", x$date.fit, "\n")
- cat("Date of summary:", x$date.summary, "\n")
-
- if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n")
-
- cat("\nEquations:\n")
- nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]])
- writeLines(strwrap(nice_diffs, exdent = 11))
- df <- x$df
- rdf <- df[2]
-
- cat("\nModel predictions using solution type", x$solution_type, "\n")
-
- cat("\nFitted using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n")
-
- if (!is.null(x$err_mod)) {
- cat("\nError model: ")
- cat(switch(x$err_mod,
- const = "Constant variance",
- obs = "Variance unique to each observed variable",
- tc = "Two-component variance function"), "\n")
-
- cat("\nError model algorithm:", x$error_model_algorithm, "\n")
- if (!is.null(x$d_3_message)) cat(x$d_3_message, "\n")
- }
-
- cat("\nStarting values for parameters to be optimised:\n")
- print(x$start)
-
- cat("\nStarting values for the transformed parameters actually optimised:\n")
- print(x$start_transformed)
-
- cat("\nFixed parameter values:\n")
- if(length(x$fixed$value) == 0) cat("None\n")
- else print(x$fixed)
-
- cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
- print(signif(x$par, digits = digits))
-
- if (x$calls > 0) {
- cat("\nParameter correlation:\n")
- if (!is.null(x$covar)){
- print(x$Corr, digits = digits, ...)
- } else {
- cat("No covariance matrix")
- }
- }
-
- cat("\nBacktransformed parameters:\n")
- cat("Confidence intervals for internally transformed parameters are asymmetric.\n")
- if ((x$version) < "0.9-36") {
- cat("To get the usual (questionable) t-test, upgrade mkin and repeat the fit.\n")
- print(signif(x$bpar, digits = digits))
- } else {
- cat("t-test (unrealistically) based on the assumption of normal distribution\n")
- cat("for estimators of untransformed parameters.\n")
- print(signif(x$bpar[, c(1, 3, 4, 5, 6)], digits = digits))
- }
-
- cat("\nFOCUS Chi2 error levels in percent:\n")
- x$errmin$err.min <- 100 * x$errmin$err.min
- print(x$errmin, digits=digits,...)
-
- printSFORB <- !is.null(x$SFORB)
- if(printSFORB){
- cat("\nEstimated Eigenvalues of SFORB model(s):\n")
- print(x$SFORB, digits=digits,...)
- }
-
- printff <- !is.null(x$ff)
- if(printff){
- cat("\nResulting formation fractions:\n")
- print(data.frame(ff = x$ff), digits=digits,...)
- }
-
- printdistimes <- !is.null(x$distimes)
- if(printdistimes){
- cat("\nEstimated disappearance times:\n")
- print(x$distimes, digits=digits,...)
- }
-
- printdata <- !is.null(x$data)
- if (printdata){
- cat("\nData:\n")
- print(format(x$data, digits = digits, ...), row.names = FALSE)
- }
-
- invisible(x)
-}
-# vim: set ts=2 sw=2 expandtab:
+if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
+
+#' Fit a kinetic model to data with one or more state variables
+#'
+#' This function maximises the likelihood of the observed data using the Port
+#' algorithm \code{\link{nlminb}}, and the specified initial or fixed
+#' parameters and starting values. In each step of the optimsation, the
+#' kinetic model is solved using the function \code{\link{mkinpredict}}. The
+#' parameters of the selected error model are fitted simultaneously with the
+#' degradation model parameters, as both of them are arguments of the
+#' likelihood function.
+#'
+#' Per default, parameters in the kinetic models are internally transformed in
+#' order to better satisfy the assumption of a normal distribution of their
+#' estimators.
+#'
+#' @param mkinmod A list of class \code{\link{mkinmod}}, containing the kinetic
+#' model to be fitted to the data, or one of the shorthand names ("SFO",
+#' "FOMC", "DFOP", "HS", "SFORB", "IORE"). If a shorthand name is given, a
+#' parent only degradation model is generated for the variable with the
+#' highest value in \code{observed}.
+#' @param observed A dataframe with the observed data. The first column called
+#' "name" must contain the name of the observed variable for each data point.
+#' The second column must contain the times of observation, named "time".
+#' The third column must be named "value" and contain the observed values.
+#' Zero values in the "value" column will be removed, with a warning, in
+#' order to avoid problems with fitting the two-component error model. This
+#' is not expected to be a problem, because in general, values of zero are
+#' not observed in degradation data, because there is a lower limit of
+#' detection.
+#' @param parms.ini A named vector of initial values for the parameters,
+#' including parameters to be optimised and potentially also fixed parameters
+#' as indicated by \code{fixed_parms}. If set to "auto", initial values for
+#' rate constants are set to default values. Using parameter names that are
+#' not in the model gives an error.
+#'
+#' It is possible to only specify a subset of the parameters that the model
+#' needs. You can use the parameter lists "bparms.ode" from a previously
+#' fitted model, which contains the differential equation parameters from
+#' this model. This works nicely if the models are nested. An example is
+#' given below.
+#' @param state.ini A named vector of initial values for the state variables of
+#' the model. In case the observed variables are represented by more than one
+#' model variable, the names will differ from the names of the observed
+#' variables (see \code{map} component of \code{\link{mkinmod}}). The default
+#' is to set the initial value of the first model variable to the mean of the
+#' time zero values for the variable with the maximum observed value, and all
+#' others to 0. If this variable has no time zero observations, its initial
+#' value is set to 100.
+#' @param err.ini A named vector of initial values for the error model
+#' parameters to be optimised. If set to "auto", initial values are set to
+#' default values. Otherwise, inital values for all error model parameters
+#' must be given.
+#' @param fixed_parms The names of parameters that should not be optimised but
+#' rather kept at the values specified in \code{parms.ini}.
+#' @param fixed_initials The names of model variables for which the initial
+#' state at time 0 should be excluded from the optimisation. Defaults to all
+#' state variables except for the first one.
+#' @param from_max_mean If this is set to TRUE, and the model has only one
+#' observed variable, then data before the time of the maximum observed value
+#' (after averaging for each sampling time) are discarded, and this time is
+#' subtracted from all remaining time values, so the time of the maximum
+#' observed mean value is the new time zero.
+#' @param solution_type If set to "eigen", the solution of the system of
+#' differential equations is based on the spectral decomposition of the
+#' coefficient matrix in cases that this is possible. If set to "deSolve", a
+#' numerical ode solver from package \code{\link{deSolve}} is used. If set to
+#' "analytical", an analytical solution of the model is used. This is only
+#' implemented for simple degradation experiments with only one state
+#' variable, i.e. with no metabolites. The default is "auto", which uses
+#' "analytical" if possible, otherwise "deSolve" if a compiler is present,
+#' and "eigen" if no compiler is present and the model can be expressed using
+#' eigenvalues and eigenvectors. This argument is passed on to the helper
+#' function \code{\link{mkinpredict}}.
+#' @param method.ode The solution method passed via \code{\link{mkinpredict}}
+#' to \code{\link{ode}} in case the solution type is "deSolve". The default
+#' "lsoda" is performant, but sometimes fails to converge.
+#' @param use_compiled If set to \code{FALSE}, no compiled version of the
+#' \code{\link{mkinmod}} model is used in the calls to
+#' \code{\link{mkinpredict}} even if a compiled version is present.
+#' @param control A list of control arguments passed to \code{\link{nlminb}}.
+#' @param transform_rates Boolean specifying if kinetic rate constants should
+#' be transformed in the model specification used in the fitting for better
+#' compliance with the assumption of normal distribution of the estimator. If
+#' TRUE, also alpha and beta parameters of the FOMC model are
+#' log-transformed, as well as k1 and k2 rate constants for the DFOP and HS
+#' models and the break point tb of the HS model. If FALSE, zero is used as
+#' a lower bound for the rates in the optimisation.
+#' @param transform_fractions Boolean specifying if formation fractions
+#' constants should be transformed in the model specification used in the
+#' fitting for better compliance with the assumption of normal distribution
+#' of the estimator. The default (TRUE) is to do transformations. If TRUE,
+#' the g parameter of the DFOP and HS models are also transformed, as they
+#' can also be seen as compositional data. The transformation used for these
+#' transformations is the \code{\link{ilr}} transformation.
+#' @param quiet Suppress printing out the current value of the negative
+#' log-likelihood after each improvement?
+#' @param atol Absolute error tolerance, passed to \code{\link{ode}}. Default
+#' is 1e-8, lower than in \code{\link{lsoda}}.
+#' @param rtol Absolute error tolerance, passed to \code{\link{ode}}. Default
+#' is 1e-10, much lower than in \code{\link{lsoda}}.
+#' @param n.outtimes The length of the dataseries that is produced by the model
+#' prediction function \code{\link{mkinpredict}}. This impacts the accuracy
+#' of the numerical solver if that is used (see \code{solution_type}
+#' argument. The default value is 100.
+#' @param error_model If the error model is "const", a constant standard
+#' deviation is assumed.
+#'
+#' If the error model is "obs", each observed variable is assumed to have its
+#' own variance.
+#'
+#' If the error model is "tc" (two-component error model), a two component
+#' error model similar to the one described by Rocke and Lorenzato (1995) is
+#' used for setting up the likelihood function. Note that this model
+#' deviates from the model by Rocke and Lorenzato, as their model implies
+#' that the errors follow a lognormal distribution for large values, not a
+#' normal distribution as assumed by this method.
+#' @param error_model_algorithm If "auto", the selected algorithm depends on
+#' the error model. If the error model is "const", unweighted nonlinear
+#' least squares fitting ("OLS") is selected. If the error model is "obs", or
+#' "tc", the "d_3" algorithm is selected.
+#'
+#' The algorithm "d_3" will directly minimize the negative log-likelihood and
+#' - independently - also use the three step algorithm described below. The
+#' fit with the higher likelihood is returned.
+#'
+#' The algorithm "direct" will directly minimize the negative log-likelihood.
+#'
+#' The algorithm "twostep" will minimize the negative log-likelihood after an
+#' initial unweighted least squares optimisation step.
+#'
+#' The algorithm "threestep" starts with unweighted least squares, then
+#' optimizes only the error model using the degradation model parameters
+#' found, and then minimizes the negative log-likelihood with free
+#' degradation and error model parameters.
+#'
+#' The algorithm "fourstep" starts with unweighted least squares, then
+#' optimizes only the error model using the degradation model parameters
+#' found, then optimizes the degradation model again with fixed error model
+#' parameters, and finally minimizes the negative log-likelihood with free
+#' degradation and error model parameters.
+#'
+#' The algorithm "IRLS" (Iteratively Reweighted Least Squares) starts with
+#' unweighted least squares, and then iterates optimization of the error
+#' model parameters and subsequent optimization of the degradation model
+#' using those error model parameters, until the error model parameters
+#' converge.
+#' @param reweight.tol Tolerance for the convergence criterion calculated from
+#' the error model parameters in IRLS fits.
+#' @param reweight.max.iter Maximum number of iterations in IRLS fits.
+#' @param trace_parms Should a trace of the parameter values be listed?
+#' @param \dots Further arguments that will be passed on to
+#' \code{\link{deSolve}}.
+#' @importFrom stats nlminb aggregate dist
+#' @return A list with "mkinfit" in the class attribute. A summary can be
+#' obtained by \code{\link{summary.mkinfit}}.
+#' @note When using the "IORE" submodel for metabolites, fitting with
+#' "transform_rates = TRUE" (the default) often leads to failures of the
+#' numerical ODE solver. In this situation it may help to switch off the
+#' internal rate transformation.
+#' @author Johannes Ranke
+#' @seealso Plotting methods \code{\link{plot.mkinfit}} and
+#' \code{\link{mkinparplot}}.
+#'
+#' Comparisons of models fitted to the same data can be made using
+#' \code{\link{AIC}} by virtue of the method \code{\link{logLik.mkinfit}}.
+#'
+#' Fitting of several models to several datasets in a single call to
+#' \code{\link{mmkin}}.
+#' @source Rocke, David M. und Lorenzato, Stefan (1995) A two-component model
+#' for measurement error in analytical chemistry. Technometrics 37(2), 176-184.
+#' @examples
+#'
+#' # Use shorthand notation for parent only degradation
+#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
+#' summary(fit)
+#'
+#' # One parent compound, one metabolite, both single first order.
+#' # Use mkinsub for convenience in model formulation. Pathway to sink included per default.
+#' SFO_SFO <- mkinmod(
+#' parent = mkinsub("SFO", "m1"),
+#' m1 = mkinsub("SFO"))
+#' # Fit the model to the FOCUS example dataset D using defaults
+#' print(system.time(fit <- mkinfit(SFO_SFO, FOCUS_2006_D,
+#' solution_type = "eigen", quiet = TRUE)))
+#' coef(fit)
+#' endpoints(fit)
+#' \dontrun{
+#' # deSolve is slower when no C compiler (gcc) was available during model generation
+#' print(system.time(fit.deSolve <- mkinfit(SFO_SFO, FOCUS_2006_D,
+#' solution_type = "deSolve")))
+#' coef(fit.deSolve)
+#' endpoints(fit.deSolve)
+#' }
+#'
+#' # Use stepwise fitting, using optimised parameters from parent only fit, FOMC
+#' \dontrun{
+#' FOMC_SFO <- mkinmod(
+#' parent = mkinsub("FOMC", "m1"),
+#' m1 = mkinsub("SFO"))
+#' # Fit the model to the FOCUS example dataset D using defaults
+#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE)
+#' # Use starting parameters from parent only FOMC fit
+#' fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, quiet = TRUE)
+#' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE,
+#' parms.ini = fit.FOMC$bparms.ode)
+#'
+#' # Use stepwise fitting, using optimised parameters from parent only fit, SFORB
+#' SFORB_SFO <- mkinmod(
+#' parent = list(type = "SFORB", to = "m1", sink = TRUE),
+#' m1 = list(type = "SFO"))
+#' # Fit the model to the FOCUS example dataset D using defaults
+#' fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, quiet = TRUE)
+#' fit.SFORB_SFO.deSolve <- mkinfit(SFORB_SFO, FOCUS_2006_D, solution_type = "deSolve",
+#' quiet = TRUE)
+#' # Use starting parameters from parent only SFORB fit (not really needed in this case)
+#' fit.SFORB = mkinfit("SFORB", FOCUS_2006_D, quiet = TRUE)
+#' fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode, quiet = TRUE)
+#' }
+#'
+#' \dontrun{
+#' # Weighted fits, including IRLS
+#' SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"),
+#' m1 = mkinsub("SFO"), use_of_ff = "max")
+#' f.noweight <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE)
+#' summary(f.noweight)
+#' f.obs <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "obs", quiet = TRUE)
+#' summary(f.obs)
+#' f.tc <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, error_model = "tc", quiet = TRUE)
+#' summary(f.tc)
+#' }
+#'
+#'
+#' @export
+mkinfit <- function(mkinmod, observed,
+ parms.ini = "auto",
+ state.ini = "auto",
+ err.ini = "auto",
+ fixed_parms = NULL,
+ fixed_initials = names(mkinmod$diffs)[-1],
+ from_max_mean = FALSE,
+ solution_type = c("auto", "analytical", "eigen", "deSolve"),
+ method.ode = "lsoda",
+ use_compiled = "auto",
+ control = list(eval.max = 300, iter.max = 200),
+ transform_rates = TRUE,
+ transform_fractions = TRUE,
+ quiet = FALSE,
+ atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
+ error_model = c("const", "obs", "tc"),
+ error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", "OLS"),
+ reweight.tol = 1e-8, reweight.max.iter = 10,
+ trace_parms = FALSE,
+ ...)
+{
+ # Check mkinmod and generate a model for the variable whith the highest value
+ # if a suitable string is given
+ parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic")
+ if (class(mkinmod) != "mkinmod") {
+ presumed_parent_name = observed[which.max(observed$value), "name"]
+ if (mkinmod[[1]] %in% parent_models_available) {
+ speclist <- list(list(type = mkinmod, sink = TRUE))
+ names(speclist) <- presumed_parent_name
+ mkinmod <- mkinmod(speclist = speclist)
+ } else {
+ stop("Argument mkinmod must be of class mkinmod or a string containing one of\n ",
+ paste(parent_models_available, collapse = ", "))
+ }
+ }
+
+ # Get the names of the state variables in the model
+ mod_vars <- names(mkinmod$diffs)
+
+ # Get the names of observed variables
+ obs_vars <- names(mkinmod$spec)
+
+ # Subset observed data with names of observed data in the model and remove NA values
+ observed <- subset(observed, name %in% obs_vars)
+ observed <- subset(observed, !is.na(value))
+
+ # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model)
+ if (any(observed$value == 0)) {
+ warning("Observations with value of zero were removed from the data")
+ observed <- subset(observed, value != 0)
+ }
+
+ # Obtain data for decline from maximum mean value if requested
+ if (from_max_mean) {
+ # This is only used for simple decline models
+ if (length(obs_vars) > 1)
+ stop("Decline from maximum is only implemented for models with a single observed variable")
+
+ means <- aggregate(value ~ time, data = observed, mean, na.rm=TRUE)
+ t_of_max <- means[which.max(means$value), "time"]
+ observed <- subset(observed, time >= t_of_max)
+ observed$time <- observed$time - t_of_max
+ }
+
+ # Number observations used for fitting
+ n_observed <- nrow(observed)
+
+ # Define starting values for parameters where not specified by the user
+ if (parms.ini[[1]] == "auto") parms.ini = vector()
+
+ # Warn for inital parameter specifications that are not in the model
+ wrongpar.names <- setdiff(names(parms.ini), mkinmod$parms)
+ if (length(wrongpar.names) > 0) {
+ warning("Initial parameter(s) ", paste(wrongpar.names, collapse = ", "),
+ " not used in the model")
+ parms.ini <- parms.ini[setdiff(names(parms.ini), wrongpar.names)]
+ }
+
+ # Warn that the sum of formation fractions may exceed one if they are not
+ # fitted in the transformed way
+ if (mkinmod$use_of_ff == "max" & transform_fractions == FALSE) {
+ warning("The sum of formation fractions may exceed one if you do not use ",
+ "transform_fractions = TRUE." )
+ for (box in mod_vars) {
+ # Stop if formation fractions are not transformed and we have no sink
+ if (mkinmod$spec[[box]]$sink == FALSE) {
+ stop("If formation fractions are not transformed during the fitting, ",
+ "it is not supported to turn off pathways to sink.\n ",
+ "Consider turning on the transformation of formation fractions or ",
+ "setting up a model with use_of_ff = 'min'.\n")
+ }
+ }
+ }
+
+ # Do not allow fixing formation fractions if we are using the ilr transformation,
+ # this is not supported
+ if (transform_fractions == TRUE && length(fixed_parms) > 0) {
+ if (any(grepl("^f_", fixed_parms))) {
+ stop("Fixing formation fractions is not supported when using the ilr ",
+ "transformation.")
+ }
+ }
+
+ # Set initial parameter values, including a small increment (salt)
+ # to avoid linear dependencies (singular matrix) in Eigenvalue based solutions
+ k_salt = 0
+ defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini))
+ for (parmname in defaultpar.names) {
+ # Default values for rate constants, depending on the parameterisation
+ if (grepl("^k", parmname)) {
+ parms.ini[parmname] = 0.1 + k_salt
+ k_salt = k_salt + 1e-4
+ }
+ # Default values for rate constants for reversible binding
+ if (grepl("free_bound$", parmname)) parms.ini[parmname] = 0.1
+ if (grepl("bound_free$", parmname)) parms.ini[parmname] = 0.02
+ # Default values for IORE exponents
+ if (grepl("^N", parmname)) parms.ini[parmname] = 1.1
+ # Default values for the FOMC, DFOP and HS models
+ if (parmname == "alpha") parms.ini[parmname] = 1
+ if (parmname == "beta") parms.ini[parmname] = 10
+ if (parmname == "k1") parms.ini[parmname] = 0.1
+ if (parmname == "k2") parms.ini[parmname] = 0.01
+ if (parmname == "tb") parms.ini[parmname] = 5
+ if (parmname == "g") parms.ini[parmname] = 0.5
+ if (parmname == "kmax") parms.ini[parmname] = 0.1
+ if (parmname == "k0") parms.ini[parmname] = 0.0001
+ if (parmname == "r") parms.ini[parmname] = 0.2
+ }
+ # Default values for formation fractions in case they are present
+ for (box in mod_vars) {
+ f_names <- mkinmod$parms[grep(paste0("^f_", box), mkinmod$parms)]
+ if (length(f_names) > 0) {
+ # We need to differentiate between default and specified fractions
+ # and set the unspecified to 1 - sum(specified)/n_unspecified
+ f_default_names <- intersect(f_names, defaultpar.names)
+ f_specified_names <- setdiff(f_names, defaultpar.names)
+ sum_f_specified = sum(parms.ini[f_specified_names])
+ if (sum_f_specified > 1) {
+ stop("Starting values for the formation fractions originating from ",
+ box, " sum up to more than 1.")
+ }
+ if (mkinmod$spec[[box]]$sink) n_unspecified = length(f_default_names) + 1
+ else {
+ n_unspecified = length(f_default_names)
+ }
+ parms.ini[f_default_names] <- (1 - sum_f_specified) / n_unspecified
+ }
+ }
+
+ # Set default for state.ini if appropriate
+ parent_name = names(mkinmod$spec)[[1]]
+ if (state.ini[1] == "auto") {
+ parent_time_0 = subset(observed, time == 0 & name == parent_name)$value
+ parent_time_0_mean = mean(parent_time_0, na.rm = TRUE)
+ if (is.na(parent_time_0_mean)) {
+ state.ini = c(100, rep(0, length(mkinmod$diffs) - 1))
+ } else {
+ state.ini = c(parent_time_0_mean, rep(0, length(mkinmod$diffs) - 1))
+ }
+ }
+
+ # Name the inital state variable values if they are not named yet
+ if(is.null(names(state.ini))) names(state.ini) <- mod_vars
+
+ # Transform initial parameter values for fitting
+ transparms.ini <- transform_odeparms(parms.ini, mkinmod,
+ transform_rates = transform_rates,
+ transform_fractions = transform_fractions)
+
+ # Parameters to be optimised:
+ # Kinetic parameters in parms.ini whose names are not in fixed_parms
+ parms.fixed <- parms.ini[fixed_parms]
+ parms.optim <- parms.ini[setdiff(names(parms.ini), fixed_parms)]
+
+ transparms.fixed <- transform_odeparms(parms.fixed, mkinmod,
+ transform_rates = transform_rates,
+ transform_fractions = transform_fractions)
+ transparms.optim <- transform_odeparms(parms.optim, mkinmod,
+ transform_rates = transform_rates,
+ transform_fractions = transform_fractions)
+
+ # Inital state variables in state.ini whose names are not in fixed_initials
+ state.ini.fixed <- state.ini[fixed_initials]
+ state.ini.optim <- state.ini[setdiff(names(state.ini), fixed_initials)]
+
+ # Preserve names of state variables before renaming initial state variable
+ # parameters
+ state.ini.optim.boxnames <- names(state.ini.optim)
+ state.ini.fixed.boxnames <- names(state.ini.fixed)
+ if(length(state.ini.optim) > 0) {
+ names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
+ }
+ if(length(state.ini.fixed) > 0) {
+ names(state.ini.fixed) <- paste(names(state.ini.fixed), "0", sep="_")
+ }
+
+ # Decide if the solution of the model can be based on a simple analytical
+ # formula, the spectral decomposition of the matrix (fundamental system)
+ # or a numeric ode solver from the deSolve package
+ # Prefer deSolve over eigen if a compiled model is present and use_compiled
+ # is not set to FALSE
+ solution_type = match.arg(solution_type)
+ if (solution_type == "analytical" && length(mkinmod$spec) > 1)
+ stop("Analytical solution not implemented for models with metabolites.")
+ if (solution_type == "eigen" && !is.matrix(mkinmod$coefmat))
+ stop("Eigenvalue based solution not possible, coefficient matrix not present.")
+ if (solution_type == "auto") {
+ if (length(mkinmod$spec) == 1) {
+ solution_type = "analytical"
+ } else {
+ if (!is.null(mkinmod$cf) & use_compiled[1] != FALSE) {
+ solution_type = "deSolve"
+ } else {
+ if (is.matrix(mkinmod$coefmat)) {
+ solution_type = "eigen"
+ if (max(observed$value, na.rm = TRUE) < 0.1) {
+ stop("The combination of small observed values (all < 0.1) and solution_type = eigen is error-prone")
+ }
+ } else {
+ solution_type = "deSolve"
+ }
+ }
+ }
+ }
+
+ # Get the error model and the algorithm for fitting
+ err_mod <- match.arg(error_model)
+ error_model_algorithm = match.arg(error_model_algorithm)
+ if (error_model_algorithm == "OLS") {
+ if (err_mod != "const") stop("OLS is only appropriate for constant variance")
+ }
+ if (error_model_algorithm == "auto") {
+ error_model_algorithm = switch(err_mod,
+ const = "OLS", obs = "d_3", tc = "d_3")
+ }
+ errparm_names <- switch(err_mod,
+ "const" = "sigma",
+ "obs" = paste0("sigma_", obs_vars),
+ "tc" = c("sigma_low", "rsd_high"))
+ errparm_names_optim <- if (error_model_algorithm == "OLS") NULL else errparm_names
+
+ # Define starting values for the error model
+ if (err.ini[1] != "auto") {
+ if (!identical(names(err.ini), errparm_names)) {
+ stop("Please supply initial values for error model components ", paste(errparm_names, collapse = ", "))
+ } else {
+ errparms = err.ini
+ }
+ } else {
+ if (err_mod == "const") {
+ errparms = 3
+ }
+ if (err_mod == "obs") {
+ errparms = rep(3, length(obs_vars))
+ }
+ if (err_mod == "tc") {
+ errparms <- c(sigma_low = 0.1, rsd_high = 0.1)
+ }
+ names(errparms) <- errparm_names
+ }
+ if (error_model_algorithm == "OLS") {
+ errparms_optim <- NULL
+ } else {
+ errparms_optim <- errparms
+ }
+
+ # Define outtimes for model solution.
+ # Include time points at which observed data are available
+ outtimes = sort(unique(c(observed$time, seq(min(observed$time),
+ max(observed$time),
+ length.out = n.outtimes))))
+
+ # Define the objective function for optimisation, including (back)transformations
+ cost_function <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...)
+ {
+ assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter
+
+ # Trace parameter values if requested and if we are actually optimising
+ if(trace_parms & update_data) cat(P, "\n")
+
+ # Determine local parameter values for the cost estimation
+ if (is.numeric(fixed_degparms)) {
+ cost_degparms <- fixed_degparms
+ cost_errparms <- P
+ degparms_fixed = TRUE
+ } else {
+ degparms_fixed = FALSE
+ }
+
+ if (is.numeric(fixed_errparms)) {
+ cost_degparms <- P
+ cost_errparms <- fixed_errparms
+ errparms_fixed = TRUE
+ } else {
+ errparms_fixed = FALSE
+ }
+
+ if (OLS) {
+ cost_degparms <- P
+ cost_errparms <- numeric(0)
+ }
+
+ if (!OLS & !degparms_fixed & !errparms_fixed) {
+ cost_degparms <- P[1:(length(P) - length(errparms))]
+ cost_errparms <- P[(length(cost_degparms) + 1):length(P)]
+ }
+
+ # Initial states for t0
+ if(length(state.ini.optim) > 0) {
+ odeini <- c(cost_degparms[1:length(state.ini.optim)], state.ini.fixed)
+ names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
+ } else {
+ odeini <- state.ini.fixed
+ names(odeini) <- state.ini.fixed.boxnames
+ }
+
+ odeparms.optim <- cost_degparms[(length(state.ini.optim) + 1):length(cost_degparms)]
+
+ if (trans == TRUE) {
+ odeparms <- c(odeparms.optim, transparms.fixed)
+ parms <- backtransform_odeparms(odeparms, mkinmod,
+ transform_rates = transform_rates,
+ transform_fractions = transform_fractions)
+ } else {
+ parms <- c(odeparms.optim, parms.fixed)
+ }
+
+ # Solve the system with current parameter values
+ out <- mkinpredict(mkinmod, parms,
+ odeini, outtimes,
+ solution_type = solution_type,
+ use_compiled = use_compiled,
+ method.ode = method.ode,
+ atol = atol, rtol = rtol, ...)
+
+ out_long <- mkin_wide_to_long(out, time = "time")
+
+ if (err_mod == "const") {
+ observed$std <- if (OLS) NA else cost_errparms["sigma"]
+ }
+ if (err_mod == "obs") {
+ std_names <- paste0("sigma_", observed$name)
+ observed$std <- cost_errparms[std_names]
+ }
+ if (err_mod == "tc") {
+ tmp <- merge(observed, out_long, by = c("time", "name"))
+ tmp$name <- ordered(tmp$name, levels = obs_vars)
+ tmp <- tmp[order(tmp$name, tmp$time), ]
+ observed$std <- sqrt(cost_errparms["sigma_low"]^2 + tmp$value.y^2 * cost_errparms["rsd_high"]^2)
+ }
+
+ cost_data <- merge(observed[c("name", "time", "value", "std")], out_long,
+ by = c("name", "time"), suffixes = c(".observed", ".predicted"))
+
+ if (OLS) {
+ # Cost is the sum of squared residuals
+ cost <- with(cost_data, sum((value.observed - value.predicted)^2))
+ } else {
+ # Cost is the negative log-likelihood
+ cost <- - with(cost_data,
+ sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE)))
+ }
+
+ # We update the current cost and data during the optimisation, not
+ # during hessian calculations
+ if (update_data) {
+
+ assign("out_predicted", out_long, inherits = TRUE)
+ assign("current_data", cost_data, inherits = TRUE)
+
+ if (cost < cost.current) {
+ assign("cost.current", cost, inherits = TRUE)
+ if (!quiet) cat(ifelse(OLS, "Sum of squared residuals", "Negative log-likelihood"),
+ " at call ", calls, ": ", cost.current, "\n", sep = "")
+ }
+ }
+ return(cost)
+ }
+
+ names_optim <- c(names(state.ini.optim),
+ names(transparms.optim),
+ errparm_names_optim)
+ n_optim <- length(names_optim)
+
+ # Define lower and upper bounds other than -Inf and Inf for parameters
+ # for which no internal transformation is requested in the call to mkinfit
+ # and for optimised error model parameters
+ lower <- rep(-Inf, n_optim)
+ upper <- rep(Inf, n_optim)
+ names(lower) <- names(upper) <- names_optim
+
+ # IORE exponents are not transformed, but need a lower bound
+ index_N <- grep("^N", names(lower))
+ lower[index_N] <- 0
+
+ if (!transform_rates) {
+ index_k <- grep("^k_", names(lower))
+ lower[index_k] <- 0
+ index_k__iore <- grep("^k__iore_", names(lower))
+ lower[index_k__iore] <- 0
+ other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2", "tb", "r"), names(lower))
+ lower[other_rate_parms] <- 0
+ }
+
+ if (!transform_fractions) {
+ index_f <- grep("^f_", names(upper))
+ lower[index_f] <- 0
+ upper[index_f] <- 1
+ other_fraction_parms <- intersect(c("g"), names(upper))
+ lower[other_fraction_parms] <- 0
+ upper[other_fraction_parms] <- 1
+ }
+
+ if (err_mod == "const") {
+ if (error_model_algorithm != "OLS") {
+ lower["sigma"] <- 0
+ }
+ }
+ if (err_mod == "obs") {
+ index_sigma <- grep("^sigma_", names(lower))
+ lower[index_sigma] <- 0
+ }
+ if (err_mod == "tc") {
+ lower["sigma_low"] <- 0
+ lower["rsd_high"] <- 0
+ }
+
+ # Counter for cost function evaluations
+ calls = 0
+ cost.current <- Inf
+ out_predicted <- NA
+ current_data <- NA
+
+ # Show parameter names if tracing is requested
+ if(trace_parms) cat(names_optim, "\n")
+
+ # browser()
+
+ # Do the fit and take the time until the hessians are calculated
+ fit_time <- system.time({
+ degparms <- c(state.ini.optim, transparms.optim)
+ n_degparms <- length(degparms)
+ degparms_index <- seq(1, n_degparms)
+ errparms_index <- seq(n_degparms + 1, length.out = length(errparms))
+
+ if (error_model_algorithm == "d_3") {
+ if (!quiet) message("Directly optimising the complete model")
+ parms.start <- c(degparms, errparms)
+ fit_direct <- nlminb(parms.start, cost_function,
+ lower = lower[names(parms.start)],
+ upper = upper[names(parms.start)],
+ control = control, ...)
+ fit_direct$logLik <- - cost.current
+ if (error_model_algorithm == "direct") {
+ degparms <- fit_direct$par[degparms_index]
+ errparms <- fit_direct$par[errparms_index]
+ } else {
+ cost.current <- Inf # reset to avoid conflict with the OLS step
+ }
+ }
+ if (error_model_algorithm != "direct") {
+ if (!quiet) message("Ordinary least squares optimisation")
+ fit <- nlminb(degparms, cost_function, control = control,
+ lower = lower[names(degparms)],
+ upper = upper[names(degparms)], OLS = TRUE, ...)
+ degparms <- fit$par
+
+ # Get the maximum likelihood estimate for sigma at the optimum parameter values
+ current_data$residual <- current_data$value.observed - current_data$value.predicted
+ sigma_mle <- sqrt(sum(current_data$residual^2)/nrow(current_data))
+
+ # Use that estimate for the constant variance, or as first guess if err_mod = "obs"
+ if (err_mod != "tc") {
+ errparms[names(errparms)] <- sigma_mle
+ }
+ fit$par <- c(fit$par, errparms)
+
+ cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
+ fit$logLik <- - cost.current
+ }
+ if (error_model_algorithm %in% c("threestep", "fourstep", "d_3")) {
+ if (!quiet) message("Optimising the error model")
+ fit <- nlminb(errparms, cost_function, control = control,
+ lower = lower[names(errparms)],
+ upper = upper[names(errparms)],
+ fixed_degparms = degparms, ...)
+ errparms <- fit$par
+ }
+ if (error_model_algorithm == "fourstep") {
+ if (!quiet) message("Optimising the degradation model")
+ fit <- nlminb(degparms, cost_function, control = control,
+ lower = lower[names(degparms)],
+ upper = upper[names(degparms)],
+ fixed_errparms = errparms, ...)
+ degparms <- fit$par
+ }
+ if (error_model_algorithm %in%
+ c("direct", "twostep", "threestep", "fourstep", "d_3")) {
+ if (!quiet) message("Optimising the complete model")
+ parms.start <- c(degparms, errparms)
+ fit <- nlminb(parms.start, cost_function,
+ lower = lower[names(parms.start)],
+ upper = upper[names(parms.start)],
+ control = control, ...)
+ degparms <- fit$par[degparms_index]
+ errparms <- fit$par[errparms_index]
+ fit$logLik <- - cost.current
+
+ if (error_model_algorithm == "d_3") {
+ d_3_messages = c(
+ same = "Direct fitting and three-step fitting yield approximately the same likelihood",
+ threestep = "Three-step fitting yielded a higher likelihood than direct fitting",
+ direct = "Direct fitting yielded a higher likelihood than three-step fitting")
+ rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
+ if (rel_diff < 0.0001) {
+ if (!quiet) message(d_3_messages["same"])
+ fit$d_3_message <- d_3_messages["same"]
+ } else {
+ if (fit$logLik > fit_direct$logLik) {
+ if (!quiet) message(d_3_messages["threestep"])
+ fit$d_3_message <- d_3_messages["threestep"]
+ } else {
+ if (!quiet) message(d_3_messages["direct"])
+ fit <- fit_direct
+ fit$d_3_message <- d_3_messages["direct"]
+ }
+ }
+ }
+ }
+ if (err_mod != "const" & error_model_algorithm == "IRLS") {
+ reweight.diff <- 1
+ n.iter <- 0
+ errparms_last <- errparms
+
+ while (reweight.diff > reweight.tol &
+ n.iter < reweight.max.iter) {
+
+ if (!quiet) message("Optimising the error model")
+ fit <- nlminb(errparms, cost_function, control = control,
+ lower = lower[names(errparms)],
+ upper = upper[names(errparms)],
+ fixed_degparms = degparms, ...)
+ errparms <- fit$par
+
+ if (!quiet) message("Optimising the degradation model")
+ fit <- nlminb(degparms, cost_function, control = control,
+ lower = lower[names(degparms)],
+ upper = upper[names(degparms)],
+ fixed_errparms = errparms, ...)
+ degparms <- fit$par
+
+ reweight.diff <- dist(rbind(errparms, errparms_last))
+ errparms_last <- errparms
+
+ fit$par <- c(fit$par, errparms)
+ cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
+ fit$logLik <- - cost.current
+ }
+ }
+
+ fit$hessian <- try(numDeriv::hessian(cost_function, c(degparms, errparms), OLS = FALSE,
+ update_data = FALSE), silent = TRUE)
+
+ # Backtransform parameters
+ bparms.optim = backtransform_odeparms(fit$par, mkinmod,
+ transform_rates = transform_rates,
+ transform_fractions = transform_fractions)
+ bparms.fixed = c(state.ini.fixed, parms.fixed)
+ bparms.all = c(bparms.optim, parms.fixed)
+
+ fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.all, errparms),
+ OLS = FALSE, trans = FALSE, update_data = FALSE), silent = TRUE)
+ })
+
+ fit$error_model_algorithm <- error_model_algorithm
+
+ if (fit$convergence != 0) {
+ fit$warning = paste0("Optimisation did not converge:\n", fit$message)
+ warning(fit$warning)
+ } else {
+ if(!quiet) message("Optimisation successfully terminated.\n")
+ }
+
+ # We need to return some more data for summary and plotting
+ fit$solution_type <- solution_type
+ fit$transform_rates <- transform_rates
+ fit$transform_fractions <- transform_fractions
+ fit$reweight.tol <- reweight.tol
+ fit$reweight.max.iter <- reweight.max.iter
+ fit$control <- control
+ fit$calls <- calls
+ fit$time <- fit_time
+
+ # We also need the model for summary and plotting
+ fit$mkinmod <- mkinmod
+
+ # We need data and predictions for summary and plotting
+ fit$observed <- observed
+ fit$obs_vars <- obs_vars
+ fit$predicted <- out_predicted
+
+ # Residual sum of squares as a function of the fitted parameters
+ fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE)
+
+ # Log-likelihood with possibility to fix degparms or errparms
+ fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE) {
+ - cost_function(P, fixed_degparms = fixed_degparms,
+ fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE)
+ }
+
+ # Collect initial parameter values in three dataframes
+ fit$start <- data.frame(value = c(state.ini.optim,
+ parms.optim, errparms_optim))
+ fit$start$type = c(rep("state", length(state.ini.optim)),
+ rep("deparm", length(parms.optim)),
+ rep("error", length(errparms_optim)))
+
+ fit$start_transformed = data.frame(
+ value = c(state.ini.optim, transparms.optim, errparms_optim),
+ lower = lower,
+ upper = upper)
+
+ fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))
+ fit$fixed$type = c(rep("state", length(state.ini.fixed)),
+ rep("deparm", length(parms.fixed)))
+
+ # Sort observed, predicted and residuals
+ current_data$name <- ordered(current_data$name, levels = obs_vars)
+
+ ordered_data <- current_data[order(current_data$name, current_data$time), ]
+
+ fit$data <- data.frame(time = ordered_data$time,
+ variable = ordered_data$name,
+ observed = ordered_data$value.observed,
+ predicted = ordered_data$value.predicted)
+
+ fit$data$residual <- fit$data$observed - fit$data$predicted
+
+ fit$atol <- atol
+ fit$rtol <- rtol
+ fit$err_mod <- err_mod
+
+ # Return different sets of backtransformed parameters for summary and plotting
+ fit$bparms.optim <- bparms.optim
+ fit$bparms.fixed <- bparms.fixed
+
+ # Return ode and state parameters for further fitting
+ fit$bparms.ode <- bparms.all[mkinmod$parms]
+ fit$bparms.state <- c(bparms.all[setdiff(names(bparms.all), names(fit$bparms.ode))],
+ state.ini.fixed)
+ names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state))
+
+ fit$errparms <- errparms
+ fit$df.residual <- n_observed - length(c(degparms, errparms))
+
+ fit$date <- date()
+ fit$version <- as.character(utils::packageVersion("mkin"))
+ fit$Rversion <- paste(R.version$major, R.version$minor, sep=".")
+
+ class(fit) <- c("mkinfit", "modFit")
+ return(fit)
+}

Contact - Imprint