mkinfit. mkin

+ +

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".

+ + +

mkinfit(mkinmod, observed,
+  parms.ini = "auto",
+  state.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",
+  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. 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. +
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. +
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. +

+ +

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.44.9000 
+#> R version:       3.3.1 
+#> Date of fit:     Thu Oct  6 09:17:59 2016 
+#> Date of summary: Thu Oct  6 09:17:59 2016 
+#> 
+#> 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.153 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.252   1.176   0.914 
+#> 
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)
+## ---------------------------------------------
+
+# 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)
+## ---------------------------------------------
+
+## 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.value <- mkinfit(SFO_SFO.ff, subset(FOCUS_2006_D, value != 0), err = "value")
+# summary(f.w.value)
+## ---------------------------------------------
+
+## 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)
+## ---------------------------------------------

Author

+ + Johannes Ranke + +

+

Arguments

Value

See also

Note

Note

Examples

Author