mkinfit. mkin 0.9-41

Usage

mkinfit(mkinmod, observed,  parms.ini = "auto", state.ini = "auto", fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1],  solution_type = c("auto", "analytical", "eigen", "deSolve"), method.ode = "lsoda", use_compiled = "auto", method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"), maxit.modFit = "auto", control.modFit = list(), transform_rates = TRUE, transform_fractions = TRUE, plot = FALSE, quiet = FALSE, err = NULL, weight = "none",  scaleVar = FALSE,  atol = 1e-8, rtol = 1e-10, n.outtimes = 100,  reweight.method = NULL, reweight.tol = 1e-8, reweight.max.iter = 10, trace_parms = FALSE, ...)

+ +

Arguments

mkinmod: + A list of class mkinmod, containing the kinetic model to be + fitted to the data, or one of the shorthand names ("SFO", "FOMC", "DFOP", + "HS", "SFORB"). If a shorthand name is given, a parent only degradation + model is generated for the variable with the highest value in + observed. +
observed: + The observed data. It has to be in the long format as described in + modFit, i.e. 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. Optionally, a further column + can contain weights for each data point. If it is not named "err", its name + must be passed as a further argument named err which is then passed + on to modFit. +
parms.ini: + A named vector of initial values for the parameters, including parameters + to be optimised and potentially also fixed parameters as indicated by + 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. +
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 map component of 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. +
fixed_parms: + The names of parameters that should not be optimised but rather kept at the + values specified in parms.ini. +
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. +
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 + 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 "eigen" if the model can be expressed using eigenvalues and + eigenvectors, and finally "deSolve" for the remaining models (time + dependence of degradation rates and metabolites). This argument is passed + on to the helper function mkinpredict. +
method.ode: + The solution method passed via mkinpredict to + ode in case the solution type is "deSolve". The default + "lsoda" is performant, but sometimes fails to converge. +
use_compiled: + If set to FALSE, no compiled version of the mkinmod + model is used, in the calls to mkinpredict even if + a compiled verion is present. +
method.modFit: + The optimisation method passed to modFit. + + In order to optimally deal with problems where local minima occur, the + "Port" algorithm is now used per default as it is less prone to get trapped + in local minima and depends less on starting values for parameters than + the Levenberg Marquardt variant selected by "Marq". However, "Port" needs + more iterations. + + The former default "Marq" is the Levenberg Marquardt algorithm + nls.lm from the package minpack.lm and usually needs + the least number of iterations. + + The "Pseudo" algorithm is not included because it needs finite parameter bounds + which are currently not supported. + + The "Newton" algorithm is not included because its number of iterations + can not be controlled by control.modFit and it does not appear + to provide advantages over the other algorithms. +
maxit.modFit: + Maximum number of iterations in the optimisation. If not "auto", this will + be passed to the method called by modFit, overriding + what may be specified in the next argument control.modFit. +
control.modFit: + Additional arguments passed to the optimisation method used by + modFit. +
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. +
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 + ilr transformation. +
plot: + Should the observed values and the numerical solutions be plotted at each + stage of the optimisation? +
quiet: + Suppress printing out the current model cost after each improvement? +
err: either NULL, or the name of the column with the + error estimates, used to weigh the residuals (see details of + modCost); if NULL, then the residuals are not weighed. +
weight: + only if err=NULL: how to weight the residuals, one of "none", + "std", "mean", see details of modCost. +
scaleVar: + Will be passed to modCost. Default is not to scale Variables + according to the number of observations. +
atol: + Absolute error tolerance, passed to ode. Default is 1e-8, + lower than in lsoda. +
rtol: + Absolute error tolerance, passed to ode. Default is 1e-10, + much lower than in lsoda. +
n.outtimes: + The length of the dataseries that is produced by the model prediction + function mkinpredict. This impacts the accuracy of + the numerical solver if that is used (see solution_type argument. + The default value is 100. +
reweight.method: + The method used for iteratively reweighting residuals, also known + as iteratively reweighted least squares (IRLS). Default is NULL, + the other method implemented is called "obs", meaning that each + observed variable is assumed to have its own variance, this is + estimated from the fit and used for weighting the residuals + in each iteration until convergence of this estimate up to + reweight.tol or up to the maximum number of iterations + specified by reweight.max.iter. +
reweight.tol: + Tolerance for convergence criterion for the variance components + in IRLS fits. +
reweight.max.iter: + Maximum iterations in IRLS fits. +
trace_parms: + Should a trace of the parameter values be listed? +
...: + Further arguments that will be passed to modFit. +

+ +

Description

+ +

This function uses the Flexible Modelling Environment package + FME to create a function calculating the model cost, i.e. the + deviation between the kinetic model and the observed data. This model cost is + then minimised using the Port algorithm nlminb, + using the specified initial or fixed parameters and starting values. + Per default, parameters in the kinetic models are internally transformed in order + to better satisfy the assumption of a normal distribution of their estimators. + In each step of the optimsation, the kinetic model is solved using the + function mkinpredict. The variance of the residuals for each + observed variable can optionally be iteratively reweighted until convergence + using the argument reweight.method = "obs".

+ +

Value

+ +

+ A list with "mkinfit" and "modFit" in the class attribute. + A summary can be obtained by summary.mkinfit. +

+ +

Note

+ +

The implementation of iteratively reweighted least squares is inspired by the + work of the KinGUII team at Bayer Crop Science (Walter Schmitt and Zhenglei + Gao). A similar implemention can also be found in CAKE 2.0, which is the + other GUI derivative of mkin, sponsored by Syngenta.

+ +

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.

+ +

Examples

# Use shorthand notation for parent only degradation
+fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
+summary(fit)
+
+mkin version:    0.9.41 
+R version:       3.2.2 
+Date of fit:     Mon Nov  9 08:25:21 2015 
+Date of summary: Mon Nov  9 08:25:21 2015 
+
+Equations:
+d_parent = - (alpha/beta) * 1/((time/beta) + 1) * parent
+
+Model predictions using solution type analytical 
+
+Fitted with method Port using 64 model solutions performed in 0.18 s
+
+Weighting: none
+
+Starting values for parameters to be optimised:
+         value   type
+parent_0  85.1  state
+alpha      1.0 deparm
+beta      10.0 deparm
+
+Starting values for the transformed parameters actually optimised:
+              value lower upper
+parent_0  85.100000  -Inf   Inf
+log_alpha  0.000000  -Inf   Inf
+log_beta   2.302585  -Inf   Inf
+
+Fixed parameter values:
+None
+
+Optimised, transformed parameters with symmetric confidence intervals:
+          Estimate Std. Error    Lower   Upper
+parent_0  85.87000     2.2460 80.38000 91.3700
+log_alpha  0.05192     0.1605 -0.34080  0.4446
+log_beta   0.65100     0.2801 -0.03452  1.3360
+
+Parameter correlation:
+          parent_0 log_alpha log_beta
+parent_0    1.0000   -0.2033  -0.3624
+log_alpha  -0.2033    1.0000   0.9547
+log_beta   -0.3624    0.9547   1.0000
+
+Residual standard error: 2.275 on 6 degrees of freedom
+
+Backtransformed parameters:
+Confidence intervals for internally transformed parameters are asymmetric.
+t-test (unrealistically) based on the assumption of normal distribution
+for estimators of untransformed parameters.
+         Estimate t value    Pr(>t)   Lower  Upper
+parent_0   85.870  38.230 1.069e-08 80.3800 91.370
+alpha       1.053   6.231 3.953e-04  0.7112  1.560
+beta        1.917   3.570 5.895e-03  0.9661  3.806
+
+Chi2 error levels in percent:
+         err.min n.optim df
+All data   6.657       3  6
+parent     6.657       3  6
+
+Estimated disappearance times:
+        DT50  DT90 DT50back
+parent 1.785 15.15     4.56
+
+Data:
+ time variable observed predicted residual
+    0   parent     85.1    85.875  -0.7749
+    1   parent     57.9    55.191   2.7091
+    3   parent     29.9    31.845  -1.9452
+    7   parent     14.6    17.012  -2.4124
+   14   parent      9.7     9.241   0.4590
+   28   parent      6.6     4.754   1.8460
+   63   parent      4.0     2.102   1.8977
+   91   parent      3.9     1.441   2.4590
+  119   parent      0.6     1.092  -0.4919
+
+
+# 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"))
+
+Successfully compiled differential equation model from auto-generated C code.
+# 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)))
+
+   user  system elapsed 
+  1.176   1.204   0.894 
+
+coef(fit)
+
+         parent_0 log_k_parent_sink   log_k_parent_m1     log_k_m1_sink 
+         99.59848          -3.03822          -2.98030          -5.24750 
+
+endpoints(fit)
+
+$ff
+parent_sink   parent_m1     m1_sink 
+   0.485524    0.514476    1.000000 
+
+$SFORB
+logical(0)
+
+$distimes
+             DT50      DT90
+parent   7.022929  23.32967
+m1     131.760712 437.69961
+
+
+## Not run: 
+# # 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)
+# ## End(Not run)
+
+# Use stepwise fitting, using optimised parameters from parent only fit, FOMC
+## Not run: 
+# 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)
+# # Use starting parameters from parent only FOMC fit
+# fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, plot=TRUE)
+# fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_2006_D, 
+#   parms.ini = fit.FOMC$bparms.ode, plot=TRUE)
+# 
+# # 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)
+# fit.SFORB_SFO.deSolve <- mkinfit(SFORB_SFO, FOCUS_2006_D, solution_type = "deSolve")
+# # Use starting parameters from parent only SFORB fit (not really needed in this case)
+# fit.SFORB = mkinfit("SFORB", FOCUS_2006_D)
+# fit.SFORB_SFO <- mkinfit(SFORB_SFO, FOCUS_2006_D, parms.ini = fit.SFORB$bparms.ode)
+# ## End(Not run)
+
+## Not run: 
+# # 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)
+# summary(f.noweight)
+# f.irls <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, reweight.method = "obs")
+# summary(f.irls)
+# f.w.mean <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, weight = "mean")
+# summary(f.w.mean)
+# f.w.mean.irls <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, weight = "mean",
+#                          reweight.method = "obs")
+# summary(f.w.mean.irls)
+# ## End(Not run)
+
+## Not run: 
+# # Manual weighting
+# dw <- FOCUS_2006_D
+# errors <- c(parent = 2, m1 = 1)
+# dw$err.man <- errors[FOCUS_2006_D$name]
+# f.w.man <- mkinfit(SFO_SFO.ff, dw, err = "err.man")
+# summary(f.w.man)
+# f.w.man.irls <- mkinfit(SFO_SFO.ff, dw, err = "err.man",
+#                        reweight.method = "obs")
+# summary(f.w.man.irls)
+# ## End(Not run)
+

+ + +

Author

+ + Johannes Ranke + + +

+ Fit a kinetic model to data with one or more state variables +

Usage

Arguments

Description

Value

Note

Note

Examples

See also

Author