+ + + + + -

Laboratory Data L1

+ +

Laboratory Data L1

The following code defines example dataset L1 from the FOCUS kinetics report, p. 284:

library("mkin")
+library("mkin", quietly = TRUE)
 FOCUS_2006_L1 = data.frame(
   t = rep(c(0, 1, 2, 3, 5, 7, 14, 21, 30), each = 2),
   parent = c(88.3, 91.4, 85.6, 84.5, 78.9, 77.6, 
@@ -88,20 +230,20 @@ FOCUS_2006_L1 = data.frame(
              27.7, 27.3, 10.0, 10.4, 2.9, 4.0))
 FOCUS_2006_L1_mkin <- mkin_wide_to_long(FOCUS_2006_L1)
 Here we use the assumptions of simple first order (SFO), the case of declining rate constant over time (FOMC) and the case of two different phases of the kinetics (DFOP). For a more detailed discussion of the models, please see the FOCUS kinetics report.
-Since mkin version 0.9-32 (July 2014), we can use shorthand notation like SFO for parent only degradation models. The following two lines fit the model and produce the summary report of the model fit. This covers the numerical analysis given in the FOCUS report.
-m.L1.SFO <- mkinfit("SFO", FOCUS_2006_L1_mkin, quiet=TRUE)
+Since mkin version 0.9-32 (July 2014), we can use shorthand notation like "SFO" for parent only degradation models. The following two lines fit the model and produce the summary report of the model fit. This covers the numerical analysis given in the FOCUS report.
+m.L1.SFO <- mkinfit("SFO", FOCUS_2006_L1_mkin, quiet = TRUE)
 summary(m.L1.SFO)
-## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:54 2016 
-## Date of summary: Thu Mar 24 08:30:54 2016 
+## mkin version:    0.9.43 
+## R version:       3.3.1 
+## Date of fit:     Tue Jun 28 01:22:22 2016 
+## Date of summary: Tue Jun 28 01:22:22 2016 
 ## 
 ## Equations:
 ## d_parent = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Port using 37 model solutions performed in 0.091 s
+## Fitted with method Port using 37 model solutions performed in 0.253 s
 ## 
 ## Weighting: none
 ## 
@@ -172,19 +314,22 @@ summary(m.L1.SFO)
 ##    30   parent      2.9     5.251  -2.3513
 ##    30   parent      4.0     5.251  -1.2513
 A plot of the fit is obtained with the plot function for mkinfit objects.
-plot(m.L1.SFO)
- The residual plot can be easily obtained by
+plot(m.L1.SFO, show_errmin = TRUE, main = "FOCUS L1 - SFO")
+
+The residual plot can be easily obtained by
 mkinresplot(m.L1.SFO, ylab = "Observed", xlab = "Time")
-
-For comparison, the FOMC model is fitted as well, and the chi^2 error level is checked.
+
+For comparison, the FOMC model is fitted as well, and the χ² error level is checked.
 m.L1.FOMC <- mkinfit("FOMC", FOCUS_2006_L1_mkin, quiet=TRUE)
 ## Warning in mkinfit("FOMC", FOCUS_2006_L1_mkin, quiet = TRUE): Optimisation by method Port did not converge.
 ## Convergence code is 1
+plot(m.L1.FOMC, show_errmin = TRUE, main = "FOCUS L1 - FOMC")
+
 summary(m.L1.FOMC, data = FALSE)
-## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:55 2016 
-## Date of summary: Thu Mar 24 08:30:55 2016 
+## mkin version:    0.9.43 
+## R version:       3.3.1 
+## Date of fit:     Tue Jun 28 01:22:24 2016 
+## Date of summary: Tue Jun 28 01:22:24 2016 
 ## 
 ## 
 ## Warning: Optimisation by method Port did not converge.
@@ -196,7 +341,7 @@ summary(m.L1.SFO)
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Port using 188 model solutions performed in 0.527 s
+## Fitted with method Port using 188 model solutions performed in 1.241 s
 ## 
 ## Weighting: none
 ## 
@@ -246,11 +391,12 @@ summary(m.L1.SFO)
 ## Estimated disappearance times:
 ##        DT50  DT90 DT50back
 ## parent 7.25 24.08     7.25
-Due to the higher number of parameters, and the lower number of degrees of freedom of the fit, the chi^2 error level is actually higher for the FOMC model (3.6%) than for the SFO model (3.4%). Additionally, the parameters log_alpha and log_beta internally fitted in the model have p-values for the two sided t-test of 0.18 and 0.125, and their correlation is 1.000, indicating that the model is overparameterised.
-The chi^2 error levels reported in Appendix 3 and Appendix 7 to the FOCUS kinetics report are rounded to integer percentages and partly deviate by one percentage point from the results calculated by mkin. The reason for this is not known. However, mkin gives the same chi^2 error levels as the kinfit package. Furthermore, the calculation routines of the kinfit package have been extensively compared to the results obtained by the KinGUI software, as documented in the kinfit package vignette. KinGUI is a widely used standard package in this field.
+We get a warning that the default optimisation algorithm Port did not converge, which is an indication that the model is overparameterised, i.e. contains too many parameters that are ill-defined as a consequence.
+And in fact, due to the higher number of parameters, and the lower number of degrees of freedom of the fit, the χ² error level is actually higher for the FOMC model (3.6%) than for the SFO model (3.4%). Additionally, the parameters log_alpha and log_beta internally fitted in the model have excessive confidence intervals, that span more than 25 orders of magnitude (!) when backtransformed to the scale of alpha and beta. Also, the t-test for significant difference from zero does not indicate such a significant difference, with p-values greater than 0.1, and finally, the parameter correlation of log_alpha and log_beta is 1.000, clearly indicating that the model is overparameterised.
+The χ² error levels reported in Appendix 3 and Appendix 7 to the FOCUS kinetics report are rounded to integer percentages and partly deviate by one percentage point from the results calculated by mkin. The reason for this is not known. However, mkin gives the same χ² error levels as the kinfit package and the calculation routines of the kinfit package have been extensively compared to the results obtained by the KinGUI software, as documented in the kinfit package vignette. KinGUI was the first widely used standard package in this field. Also, the calculation of χ² error levels was compared with KinGUII, CAKE and DegKin manager in a project sponsored by the German Umweltbundesamt (Ranke).

Laboratory Data L2

The following code defines example dataset L2 from the FOCUS kinetics report, p. 287:

FOCUS_2006_L2 = data.frame(
   t = rep(c(0, 1, 3, 7, 14, 28), each = 2),
@@ -258,108 +404,36 @@ summary(m.L1.SFO)

19.3, 22.3, 4.6, 4.6, 2.6, 1.2, 0.3, 0.6)) FOCUS_2006_L2_mkin <- mkin_wide_to_long(FOCUS_2006_L2) -

Again, the SFO model is fitted and a summary is obtained:

SFO fit for L2

Again, the SFO model is fitted and the result is plotted. The residual plot can be obtained simply by adding the argument show_residuals to the plot command.

m.L2.SFO <- mkinfit("SFO", FOCUS_2006_L2_mkin, quiet=TRUE)
-summary(m.L2.SFO)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:55 2016 
-## Date of summary: Thu Mar 24 08:30:55 2016 
-## 
-## Equations:
-## d_parent = - k_parent_sink * parent
-## 
-## Model predictions using solution type analytical 
-## 
-## Fitted with method Port using 41 model solutions performed in 0.099 s
-## 
-## Weighting: none
-## 
-## Starting values for parameters to be optimised:
-##               value   type
-## parent_0      93.95  state
-## k_parent_sink  0.10 deparm
-## 
-## Starting values for the transformed parameters actually optimised:
-##                       value lower upper
-## parent_0          93.950000  -Inf   Inf
-## log_k_parent_sink -2.302585  -Inf   Inf
-## 
-## Fixed parameter values:
-## None
-## 
-## Optimised, transformed parameters with symmetric confidence intervals:
-##                   Estimate Std. Error   Lower   Upper
-## parent_0           91.4700     3.8070 82.9800 99.9500
-## log_k_parent_sink  -0.4112     0.1074 -0.6505 -0.1719
-## 
-## Parameter correlation:
-##                   parent_0 log_k_parent_sink
-## parent_0            1.0000            0.4295
-## log_k_parent_sink   0.4295            1.0000
-## 
-## Residual standard error: 5.51 on 10 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       91.4700   24.03 1.773e-10 82.9800 99.9500
-## k_parent_sink   0.6629    9.31 1.525e-06  0.5218  0.8421
-## 
-## Chi2 error levels in percent:
-##          err.min n.optim df
-## All data   14.38       2  4
-## parent     14.38       2  4
-## 
-## Resulting formation fractions:
-##             ff
-## parent_sink  1
-## 
-## Estimated disappearance times:
-##         DT50  DT90
-## parent 1.046 3.474
-## 
-## Data:
-##  time variable observed predicted residual
-##     0   parent     96.1 9.147e+01   4.6343
-##     0   parent     91.8 9.147e+01   0.3343
-##     1   parent     41.4 4.714e+01  -5.7394
-##     1   parent     38.7 4.714e+01  -8.4394
-##     3   parent     19.3 1.252e+01   6.7790
-##     3   parent     22.3 1.252e+01   9.7790
-##     7   parent      4.6 8.834e-01   3.7166
-##     7   parent      4.6 8.834e-01   3.7166
-##    14   parent      2.6 8.532e-03   2.5915
-##    14   parent      1.2 8.532e-03   1.1915
-##    28   parent      0.3 7.958e-07   0.3000
-##    28   parent      0.6 7.958e-07   0.6000

The chi^2 error level of 14% suggests that the model does not fit very well. This is also obvious from the plots of the fit and the residuals.

par(mfrow = c(2, 1))
-plot(m.L2.SFO)
-mkinresplot(m.L2.SFO)

+plot(m.L2.SFO, show_residuals = TRUE, show_errmin = TRUE, + main = "FOCUS L2 - SFO") +

The χ² error level of 14% suggests that the model does not fit very well. This is also obvious from the plots of the fit, in which we have included the residual plot.

In the FOCUS kinetics report, it is stated that there is no apparent systematic error observed from the residual plot up to the measured DT90 (approximately at day 5), and there is an underestimation beyond that point.

We may add that it is difficult to judge the random nature of the residuals just from the three samplings at days 0, 1 and 3. Also, it is not clear a priori why a consistent underestimation after the approximate DT90 should be irrelevant. However, this can be rationalised by the fact that the FOCUS fate models generally only implement SFO kinetics.

For comparison, the FOMC model is fitted as well, and the chi^2 error level is checked.

FOMC fit for L2

For comparison, the FOMC model is fitted as well, and the χ² error level is checked.

m.L2.FOMC <- mkinfit("FOMC", FOCUS_2006_L2_mkin, quiet = TRUE)
-par(mfrow = c(2, 1))
-plot(m.L2.FOMC)
-mkinresplot(m.L2.FOMC)

+plot(m.L2.FOMC, show_residuals = TRUE, + main = "FOCUS L2 - FOMC") +

summary(m.L2.FOMC, data = FALSE)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:56 2016 
-## Date of summary: Thu Mar 24 08:30:56 2016 
+## mkin version:    0.9.43 
+## R version:       3.3.1 
+## Date of fit:     Tue Jun 28 01:22:27 2016 
+## Date of summary: Tue Jun 28 01:22:27 2016 
 ## 
 ## Equations:
 ## d_parent = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Port using 81 model solutions performed in 0.2 s
+## Fitted with method Port using 81 model solutions performed in 0.541 s
 ## 
 ## Weighting: none
 ## 
@@ -409,22 +483,20 @@ mkinresplot(m.L2.FOMC)
 ## Estimated disappearance times:
 ##          DT50  DT90 DT50back
 ## parent 0.8092 5.356    1.612

The error level at which the chi^2 test passes is much lower in this case. Therefore, the FOMC model provides a better description of the data, as less experimental error has to be assumed in order to explain the data.

Fitting the four parameter DFOP model further reduces the chi^2 error level.

The error level at which the χ² test passes is much lower in this case. Therefore, the FOMC model provides a better description of the data, as less experimental error has to be assumed in order to explain the data.

DFOP fit for L2

Fitting the four parameter DFOP model further reduces the χ² error level.

m.L2.DFOP <- mkinfit("DFOP", FOCUS_2006_L2_mkin, quiet = TRUE)
-plot(m.L2.DFOP)

Here, the default starting parameters for the DFOP model obviously do not lead to a reasonable solution. Therefore the fit is repeated with different starting parameters.

m.L2.DFOP <- mkinfit("DFOP", FOCUS_2006_L2_mkin, 
-  parms.ini = c(k1 = 1, k2 = 0.01, g = 0.8),
-  quiet=TRUE)
-plot(m.L2.DFOP)

+plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE, + main = "FOCUS L2 - DFOP") +

summary(m.L2.DFOP, data = FALSE)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:58 2016 
-## Date of summary: Thu Mar 24 08:30:58 2016 
+## mkin version:    0.9.43 
+## R version:       3.3.1 
+## Date of fit:     Tue Jun 28 01:22:30 2016 
+## Date of summary: Tue Jun 28 01:22:30 2016 
 ## 
 ## Equations:
 ## d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -433,23 +505,23 @@ plot(m.L2.DFOP)
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Port using 336 model solutions performed in 0.903 s
+## Fitted with method Port using 336 model solutions performed in 2.29 s
 ## 
 ## Weighting: none
 ## 
 ## Starting values for parameters to be optimised:
 ##          value   type
 ## parent_0 93.95  state
-## k1        1.00 deparm
+## k1        0.10 deparm
 ## k2        0.01 deparm
-## g         0.80 deparm
+## g         0.50 deparm
 ## 
 ## Starting values for the transformed parameters actually optimised:
-##               value lower upper
-## parent_0 93.9500000  -Inf   Inf
-## log_k1    0.0000000  -Inf   Inf
-## log_k2   -4.6051702  -Inf   Inf
-## g_ilr     0.9802581  -Inf   Inf
+##              value lower upper
+## parent_0 93.950000  -Inf   Inf
+## log_k1   -2.302585  -Inf   Inf
+## log_k2   -4.605170  -Inf   Inf
+## g_ilr     0.000000  -Inf   Inf
 ## 
 ## Fixed parameter values:
 ## None
@@ -457,12 +529,13 @@ plot(m.L2.DFOP)

## Optimised, transformed parameters with symmetric confidence intervals: ## Estimate Std. Error Lower Upper ## parent_0 93.9500 NA NA NA -## log_k1 3.1210 NA NA NA +## log_k1 3.1370 NA NA NA ## log_k2 -1.0880 NA NA NA ## g_ilr -0.2821 NA NA NA ## -## Parameter correlation: -## Could not estimate covariance matrix; singular system: +## Parameter correlation: +

## Warning in print.summary.mkinfit(x): Could not estimate covariance matrix; singular system:

## Could not estimate covariance matrix; singular system:
 ## 
 ## Residual standard error: 1.732 on 8 degrees of freedom
 ## 
@@ -472,7 +545,7 @@ plot(m.L2.DFOP)

## for estimators of untransformed parameters. ## Estimate t value Pr(>t) Lower Upper ## parent_0 93.9500 NA NA NA NA -## k1 22.6700 NA NA NA NA +## k1 23.0400 NA NA NA NA ## k2 0.3369 NA NA NA NA ## g 0.4016 NA NA NA NA ## @@ -483,167 +556,36 @@ plot(m.L2.DFOP) ## ## Estimated disappearance times: ## DT50 DT90 DT50_k1 DT50_k2 -## parent 0.5335 5.311 0.03058 2.058 +## parent 0.5335 5.311 0.03009 2.058

Here, the DFOP model is clearly the best-fit model for dataset L2 based on the chi^2 error level criterion. However, the failure to calculate the covariance matrix indicates that the parameter estimates correlate excessively. Therefore, the FOMC model may be preferred for this dataset.

Laboratory Data L3

The following code defines example dataset L3 from the FOCUS kinetics report, p. 290.

FOCUS_2006_L3 = data.frame(
   t = c(0, 3, 7, 14, 30, 60, 91, 120),
   parent = c(97.8, 60, 51, 43, 35, 22, 15, 12))
 FOCUS_2006_L3_mkin <- mkin_wide_to_long(FOCUS_2006_L3)

SFO model, summary and plot:

m.L3.SFO <- mkinfit("SFO", FOCUS_2006_L3_mkin, quiet = TRUE)
-plot(m.L3.SFO)

summary(m.L3.SFO)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:58 2016 
-## Date of summary: Thu Mar 24 08:30:58 2016 
-## 
-## Equations:
-## d_parent = - k_parent_sink * parent
-## 
-## Model predictions using solution type analytical 
-## 
-## Fitted with method Port using 43 model solutions performed in 0.104 s
-## 
-## Weighting: none
-## 
-## Starting values for parameters to be optimised:
-##               value   type
-## parent_0       97.8  state
-## k_parent_sink   0.1 deparm
-## 
-## Starting values for the transformed parameters actually optimised:
-##                       value lower upper
-## parent_0          97.800000  -Inf   Inf
-## log_k_parent_sink -2.302585  -Inf   Inf
-## 
-## Fixed parameter values:
-## None
-## 
-## Optimised, transformed parameters with symmetric confidence intervals:
-##                   Estimate Std. Error  Lower Upper
-## parent_0            74.870     8.4570 54.180 95.57
-## log_k_parent_sink   -3.678     0.3261 -4.476 -2.88
-## 
-## Parameter correlation:
-##                   parent_0 log_k_parent_sink
-## parent_0            1.0000            0.5483
-## log_k_parent_sink   0.5483            1.0000
-## 
-## Residual standard error: 12.91 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      74.87000   8.853 5.776e-05 54.18000 95.57000
-## k_parent_sink  0.02527   3.067 1.102e-02  0.01138  0.05612
-## 
-## Chi2 error levels in percent:
-##          err.min n.optim df
-## All data   21.24       2  6
-## parent     21.24       2  6
-## 
-## Resulting formation fractions:
-##             ff
-## parent_sink  1
-## 
-## Estimated disappearance times:
-##         DT50  DT90
-## parent 27.43 91.13
-## 
-## Data:
-##  time variable observed predicted residual
-##     0   parent     97.8    74.872  22.9281
-##     3   parent     60.0    69.406  -9.4061
-##     7   parent     51.0    62.734 -11.7340
-##    14   parent     43.0    52.564  -9.5638
-##    30   parent     35.0    35.084  -0.0839
-##    60   parent     22.0    16.440   5.5602
-##    91   parent     15.0     7.511   7.4887
-##   120   parent     12.0     3.610   8.3903

The chi^2 error level of 21% as well as the plot suggest that the model does not fit very well.

The FOMC model performs better:

m.L3.FOMC <- mkinfit("FOMC", FOCUS_2006_L3_mkin, quiet = TRUE)
-plot(m.L3.FOMC)

summary(m.L3.FOMC, data = FALSE)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:59 2016 
-## Date of summary: Thu Mar 24 08:30:59 2016 
-## 
-## Equations:
-## d_parent = - (alpha/beta) * 1/((time/beta) + 1) * parent
-## 
-## Model predictions using solution type analytical 
-## 
-## Fitted with method Port using 83 model solutions performed in 0.208 s
-## 
-## Weighting: none
-## 
-## Starting values for parameters to be optimised:
-##          value   type
-## parent_0  97.8  state
-## alpha      1.0 deparm
-## beta      10.0 deparm
-## 
-## Starting values for the transformed parameters actually optimised:
-##               value lower upper
-## parent_0  97.800000  -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   96.9700     4.5500 85.2800 108.7000
-## log_alpha  -0.8619     0.1704 -1.3000  -0.4238
-## log_beta    0.6193     0.4744 -0.6003   1.8390
-## 
-## Parameter correlation:
-##           parent_0 log_alpha log_beta
-## parent_0    1.0000   -0.1512  -0.4271
-## log_alpha  -0.1512    1.0000   0.9110
-## log_beta   -0.4271    0.9110   1.0000
-## 
-## Residual standard error: 4.572 on 5 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  96.9700  21.310 2.108e-06 85.2800 108.7000
-## alpha      0.4224   5.867 1.020e-03  0.2725   0.6546
-## beta       1.8580   2.108 4.444e-02  0.5487   6.2890
-## 
-## Chi2 error levels in percent:
-##          err.min n.optim df
-## All data    7.32       3  5
-## parent      7.32       3  5
-## 
-## Estimated disappearance times:
-##         DT50  DT90 DT50back
-## parent 7.729 431.2    129.8

The error level at which the chi^2 test passes is 7% in this case.

Fitting the four parameter DFOP model further reduces the chi^2 error level considerably:

m.L3.DFOP <- mkinfit("DFOP", FOCUS_2006_L3_mkin, quiet = TRUE)
-plot(m.L3.DFOP)

summary(m.L3.DFOP, data = FALSE)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:30:59 2016 
-## Date of summary: Thu Mar 24 08:30:59 2016 
+
+Use mmkin to fit multiple models
+As of mkin version 0.9-39 (June 2015), we can fit several models to one or more datasets in one call to the function mmkin. The datasets have to be passed in a list, in this case a named list holding only the L3 dataset prepared above.
+# Only use one core here, not to offend the CRAN checks
+mm.L3 <- mmkin(c("SFO", "FOMC", "DFOP"), cores = 1,
+               list("FOCUS L3" = FOCUS_2006_L3_mkin), quiet = TRUE)
+plot(mm.L3)
+
+The χ² error level of 21% as well as the plot suggest that the SFO model does not fit very well. The FOMC model performs better, with an error level at which the χ² test passes of 7%. Fitting the four parameter DFOP model further reduces the χ² error level considerably.
+
+
+Accessing elements of mmkin objects
+The objects returned by mmkin are arranged like a matrix, with models as a row index and datasets as a column index.
+We can extract the summary and plot for e.g. the DFOP fit, using square brackets for indexing which will result in the use of the summary and plot functions working on mkinfit objects.
+summary(mm.L3[["DFOP", 1]])
+## mkin version:    0.9.43 
+## R version:       3.3.1 
+## Date of fit:     Tue Jun 28 01:22:33 2016 
+## Date of summary: Tue Jun 28 01:22:33 2016 
 ## 
 ## Equations:
 ## d_parent = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -652,7 +594,7 @@ plot(m.L3.DFOP)
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Port using 137 model solutions performed in 0.351 s
+## Fitted with method Port using 137 model solutions performed in 0.959 s
 ## 
 ## Weighting: none
 ## 
@@ -706,33 +648,51 @@ plot(m.L3.DFOP)

## ## Estimated disappearance times: ## DT50 DT90 DT50_k1 DT50_k2 -## parent 7.464 123 1.343 50.37 -

Here, a look to the model plot, the confidence intervals of the parameters and the correlation matrix suggest that the parameter estimates are reliable, and the DFOP model can be used as the best-fit model based on the chi^2 error level criterion for laboratory data L3.

+## parent 7.464 123 1.343 50.37 +## +## Data: +## time variable observed predicted residual +## 0 parent 97.8 97.75 0.05396 +## 3 parent 60.0 60.45 -0.44933 +## 7 parent 51.0 49.44 1.56338 +## 14 parent 43.0 43.84 -0.83632 +## 30 parent 35.0 35.15 -0.14707 +## 60 parent 22.0 23.26 -1.25919 +## 91 parent 15.0 15.18 -0.18181 +## 120 parent 12.0 10.19 1.81395 +

plot(mm.L3[["DFOP", 1]], show_errmin = TRUE)

This is also an example where the standard t-test for the parameter g_ilr is misleading, as it tests for a significant difference from zero. In this case, zero appears to be the correct value for this parameter, and the confidence interval for the backtransformed parameter g is quite narrow.

Laboratory Data L4

The following code defines example dataset L4 from the FOCUS kinetics report, p. 293:

FOCUS_2006_L4 = data.frame(
   t = c(0, 3, 7, 14, 30, 60, 91, 120),
   parent = c(96.6, 96.3, 94.3, 88.8, 74.9, 59.9, 53.5, 49.0))
 FOCUS_2006_L4_mkin <- mkin_wide_to_long(FOCUS_2006_L4)

SFO model, summary and plot:

m.L4.SFO <- mkinfit("SFO", FOCUS_2006_L4_mkin, quiet = TRUE)
-plot(m.L4.SFO)

summary(m.L4.SFO, data = FALSE)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:31:00 2016 
-## Date of summary: Thu Mar 24 08:31:00 2016 
+Fits of the SFO and FOMC models, plots and summaries are produced below:
+# Only use one core here, not to offend the CRAN checks
+mm.L4 <- mmkin(c("SFO", "FOMC"), cores = 1,
+               list("FOCUS L4" = FOCUS_2006_L4_mkin), 
+               quiet = TRUE)
+plot(mm.L4)
+
+The χ² error level of 3.3% as well as the plot suggest that the SFO model fits very well. The error level at which the χ² test passes is slightly lower for the FOMC model. However, the difference appears negligible.
+summary(mm.L4[["SFO", 1]], data = FALSE)
+## mkin version:    0.9.43 
+## R version:       3.3.1 
+## Date of fit:     Tue Jun 28 01:22:34 2016 
+## Date of summary: Tue Jun 28 01:22:35 2016 
 ## 
 ## Equations:
 ## d_parent = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Port using 46 model solutions performed in 0.116 s
+## Fitted with method Port using 46 model solutions performed in 0.302 s
 ## 
 ## Weighting: none
 ## 
@@ -781,23 +741,18 @@ plot(m.L4.SFO)
 ## Estimated disappearance times:
 ##        DT50 DT90
 ## parent  106  352

The chi^2 error level of 3.3% as well as the plot suggest that the model fits very well.

The FOMC model for comparison:

m.L4.FOMC <- mkinfit("FOMC", FOCUS_2006_L4_mkin, quiet = TRUE)
-plot(m.L4.FOMC)

summary(m.L4.FOMC, data = FALSE)

## mkin version:    0.9.42 
-## R version:       3.2.4 
-## Date of fit:     Thu Mar 24 08:31:00 2016 
-## Date of summary: Thu Mar 24 08:31:00 2016 
+summary(mm.L4[["FOMC", 1]], data = FALSE)
+## mkin version:    0.9.43 
+## R version:       3.3.1 
+## Date of fit:     Tue Jun 28 01:22:35 2016 
+## Date of summary: Tue Jun 28 01:22:35 2016 
 ## 
 ## Equations:
 ## d_parent = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted with method Port using 66 model solutions performed in 0.164 s
+## Fitted with method Port using 66 model solutions performed in 0.438 s
 ## 
 ## Weighting: none
 ## 
@@ -847,10 +802,17 @@ plot(m.L4.FOMC)
 ## Estimated disappearance times:
 ##         DT50 DT90 DT50back
 ## parent 108.9 1644    494.9

The error level at which the chi^2 test passes is slightly lower for the FOMC model. However, the difference appears negligible.

References

Ranke, Johannes. Prüfung und Validierung von Modellierungssoftware als Alternative zu ModelMaker 4.0. Umweltbundesamt Projektnummer 27452.

+ +

+ + + + - - + + - + @@ -38,7 +42,34 @@ if (window.hljs && document.readyState && document.readyState === "complete") { - + + + @@ -54,27 +85,140 @@ code { color: inherit; background-color: rgba(0, 0, 0, 0.04); } -img { - max-width:100%; - height: auto; +img { + max-width:100%; + height: auto; +} +.tabbed-pane { + padding-top: 12px; +} +button.code-folding-btn:focus { + outline: none; } + +

+ + + -

Performance benefit by using compiled model definitions in mkin

Johannes Ranke

2016-03-24

+ + + + + + + + + + +

Benchmark for a model that can also be solved with Eigenvalues
Benchmark for a model that can not be solved with Eigenvalues

+ + + + + +

Benchmark for a model that can also be solved with Eigenvalues

This evaluation is taken from the example section of mkinfit. When using an mkin version equal to or greater than 0.9-36 and a C compiler (gcc) is available, you will see a message that the model is being compiled from autogenerated C code when defining a model using mkinmod. The mkinmod() function checks for presence of the gcc compiler using

@@ -89,10 +233,8 @@ SFO_SFO <- mkinmod(

## Successfully compiled differential equation model from auto-generated C code.

We can compare the performance of the Eigenvalue based solution against the compiled version and the R implementation of the differential equations using the microbenchmark package.

library("microbenchmark")
-library("ggplot2")

## Need help? Try the ggplot2 mailing list:
-## http://groups.google.com/group/ggplot2.

mb.1 <- microbenchmark(
+library("ggplot2")
+mb.1 <- microbenchmark(
   "deSolve, not compiled" = mkinfit(SFO_SFO, FOCUS_2006_D, 
                                     solution_type = "deSolve", 
                                     use_compiled = FALSE, quiet = TRUE),
@@ -106,23 +248,23 @@ library("ggplot2")

## precision.

smb.1 <- summary(mb.1)
 print(mb.1)

## Unit: milliseconds
+## Unit: seconds
 ##                   expr       min        lq      mean    median        uq
-##  deSolve, not compiled 9280.0854 9299.6757 9323.2559 9319.2659 9344.8411
-##       Eigenvalue based  885.7475  891.8548  907.2823  897.9621  918.0498
-##      deSolve, compiled  713.2624  721.4990  728.2856  729.7357  735.7972
+##  deSolve, not compiled 25.634645 25.717260 25.925956 25.799875 26.071612
+##       Eigenvalue based  2.235689  2.237675  2.246796  2.239661  2.252349
+##      deSolve, compiled  1.869556  1.881345  1.887495  1.893133  1.896465
 ##        max neval cld
-##  9370.4163     3   c
-##   938.1374     3  b 
-##   741.8588     3 a
+##  26.343349     3   b
+##   2.265038     3  a 
+##   1.899797     3  a

autoplot(mb.1)

We see that using the compiled model is by a factor of 12.8 faster than using the R version with the default ode solver, and it is even faster than the Eigenvalue based solution implemented in R which does not need iterative solution of the ODEs:

We see that using the compiled model is by a factor of 13.6 faster than using the R version with the default ode solver, and it is even faster than the Eigenvalue based solution implemented in R which does not need iterative solution of the ODEs:

rownames(smb.1) <- smb.1$expr
 smb.1["median"]/smb.1["deSolve, compiled", "median"]

##                          median
-## deSolve, not compiled 12.770742
-## Eigenvalue based       1.230531
+## deSolve, not compiled 13.628133
+## Eigenvalue based       1.183044
 ## deSolve, compiled      1.000000

@@ -144,26 +286,30 @@ smb.1["median"]/smb.1["deSolve, compiled", "median" print(mb.2)

## Unit: seconds
 ##                   expr       min        lq      mean    median        uq
-##  deSolve, not compiled 20.543131 20.661195 20.720383 20.779259 20.809008
-##      deSolve, compiled  1.314865  1.316439  1.328049  1.318014  1.334642
-##       max neval cld
-##  20.83876     3   b
-##   1.35127     3  a

+## deSolve, not compiled 53.993961 54.012379 54.225054 54.030797 54.340600 +## deSolve, compiled 3.450737 3.535551 3.583034 3.620364 3.649183 +## max neval cld +## 54.650403 3 b +## 3.678001 3 a

smb.2["median"]/smb.2["deSolve, compiled", "median"]

##   median
 ## 1     NA
 ## 2     NA

autoplot(mb.2)

Here we get a performance benefit of a factor of 15.8 using the version of the differential equation model compiled from C code using the inline package!

This vignette was built with mkin 0.9.42 on

## R version 3.2.4 Revised (2016-03-16 r70336)
+
+Here we get a performance benefit of a factor of 14.9 using the version of the differential equation model compiled from C code!
+This vignette was built with mkin 0.9.43 on
+## R version 3.3.1 (2016-06-21)
 ## Platform: x86_64-pc-linux-gnu (64-bit)
 ## Running under: Debian GNU/Linux 8 (jessie)
 ## CPU model: Intel(R) Core(TM) i7-4710MQ CPU @ 2.50GHz

+
+


 
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Show All
+Hide All
+
+
+
+
+mkin - Kinetic evaluation of chemical degradation data
+Johannes Ranke
+2016-06-28
+
+
+
+
+Wissenschaftlicher Berater, Kronacher Str. 8, 79639 Grenzach-Wyhlen, Germany
 Privatdozent at the University of Bremen
+
+Abstract
+In the regulatory evaluation of chemical substances like plant protection products (pesticides), biocides and other chemicals, degradation data play an important role. For the evaluation of pesticide degradation experiments, detailed guidance has been developed, based on nonlinear optimisation. The R add-on package mkin (Ranke 2016) implements fitting some of the models recommended in this guidance from within R and calculates some statistical measures for data series within one or more compartments, for parent and metabolites.
+require(mkin)
+m_SFO_SFO_SFO <- mkinmod(parent = mkinsub("SFO", "M1"),
+                         M1 = mkinsub("SFO", "M2"),
+                         M2 = mkinsub("SFO"), 
+                         use_of_ff = "max", quiet = TRUE)
+
+sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+
+d_SFO_SFO_SFO <- mkinpredict(m_SFO_SFO_SFO,
+  c(k_parent = 0.03, 
+    f_parent_to_M1 = 0.5, k_M1 = log(2)/100, 
+    f_M1_to_M2 = 0.9, k_M2 = log(2)/50),
+  c(parent = 100, M1 = 0, M2 = 0),
+  sampling_times)
+
+d_SFO_SFO_SFO_err <- add_err(d_SFO_SFO_SFO, function(x) 3, n = 1, seed = 123456789 )
+
+f_SFO_SFO_SFO <- mkinfit(m_SFO_SFO_SFO, d_SFO_SFO_SFO_err[[1]], quiet = TRUE)
+
+plot_sep(f_SFO_SFO_SFO, lpos = c("topright", "bottomright", "bottomright"))
+
+
+
+Background
+Many approaches are possible regarding the evaluation of chemical degradation data.
+The now deprecated kinfit package (Ranke 2015) in R (R Development Core Team 2016) implements the approach recommended in the kinetics report provided by the FOrum for Co-ordination of pesticide fate models and their USe (FOCUS Work Group on Degradation Kinetics 2006; 2014) for simple data series for one parent compound in one compartment.
+The mkin package (Ranke 2016) extends this approach to data series with transformation products, commonly termed metabolites, and to more than one compartment. It is also possible to include back reactions, so equilibrium reactions and equilibrium partitioning can be specified, although this oftentimes leads to an overparameterisation of the model.
+When the first mkin code was published in 2010, the most commonly used tools for fitting more complex kinetic degradation models to experimental data were KinGUI (Schäfer et al. 2007), a MATLAB based tool with a graphical user interface that was specifically tailored to the task and included some output as proposed by the FOCUS Kinetics Workgroup, and ModelMaker, a general purpose compartment based tool providing infrastructure for fitting dynamic simulation models based on differential equations to data.
+The code was first uploaded to the BerliOS platform. When this was taken down, the version control history was imported into the R-Forge site, where the code is still mirrored today (see e.g. the initial commit on 11 May 2010).
+At that time, the R package FME (Flexible Modelling Environment) (Soetaert and Petzoldt 2010) was already available, and provided a good basis for developing a package specifically tailored to the task. The remaining challenge was to make it as easy as possible for the users (including the author of this vignette) to specify the system of differential equations and to include the output requested by the FOCUS guidance, such as the relative standard deviation that has to be assumed for the residuals, such that the \(\chi^2\) goodness-of-fit test as defined by the FOCUS kinetics workgroup would pass using an significance level \(\alpha\) of 0.05.
+Also, mkin introduced using analytical solutions for parent only kinetics for improved optimization speed. Later, Eigenvalue based solutions were introduced to mkin for the case of linear differential equations (i.e. where the FOMC or DFOP models were not used for the parent compound), greatly improving the optimization speed for these cases.
+The possibility to specify back-reactions and a biphasic model (SFORB) for metabolites were present in mkin from the very beginning.
+
+Derived software tools
+Soon after the publication of mkin, two derived tools were published, namely KinGUII (available from Bayer Crop Science) and CAKE (commissioned to Tessella by Syngenta), which added a graphical user interface (GUI), and added fitting by iteratively reweighted least squares (IRLS) and characterisation of likely parameter distributions by Markov Chain Monte Carlo (MCMC) sampling.
+CAKE focuses on a smooth use experience, sacrificing some flexibility in the model definition, originally allowing only two primary metabolites in parallel. The current version 3.2 of CAKE release in March 2016 uses a basic scheme for up to six metabolites in a flexible arrangement, but does not support back-reactions (non-instantaneous equilibria) or biphasic kinetics for metabolites.
+KinGUI offers an even more flexible widget for specifying complex kinetic models. Back-reactions (non-instanteneous equilibria) were supported early on, but until 2014, only simple first-order models could be specified for transformation products. Starting with KinGUII version 2.1, biphasic modelling of metabolites was also available in KinGUII.
+A further graphical user interface (GUI) that has recently been brought to a decent degree of maturity is the browser based GUI named gmkin. Please see its documentation page and manual for further information.
+
+
+Recent developments
+Currently (June 2016), the main features available in mkin which are not present in KinGUII or CAKE, are the speed increase by using compiled code when a compiler is present, parallel model fitting on multicore machines using the mmkin function, and the estimation of parameter confidence intervals based on transformed parameters. These are explained in more detail below.
+
+
+
+Internal parameter transformations
+For rate constants, the log transformation is used, as proposed by Bates and Watts (1988, 77, 149). Approximate intervals are constructed for the transformed rate constants (compare Bates and Watts 1988, 135), i.e. for their logarithms. Confidence intervals for the rate constants are then obtained using the appropriate backtransformation using the exponential function.
+In the first version of mkin allowing for specifying models using formation fractions, a home-made reparameterisation was used in order to ensure that the sum of formation fractions would not exceed unity.
+This method is still used in the current version of KinGUII (v2.1 from April 2014), with a modification that allows for fixing the pathway to sink to zero. CAKE uses penalties in the objective function in order to enforce this constraint.
+In 2012, an alternative reparameterisation of the formation fractions was proposed together with René Lehmann (Ranke and Lehmann 2012), based on isometric logratio transformation (ILR). The aim was to improve the validity of the linear approximation of the objective function during the parameter estimation procedure as well as in the subsequent calculation of parameter confidence intervals.
+
+Confidence intervals based on transformed parameters
+In the first attempt at providing improved parameter confidence intervals introduced to mkin in 2013, confidence intervals obtained from FME on the transformed parameters were simply all backtransformed one by one to yield asymetric confidence intervals for the backtransformed parameters.
+However, while there is a 1:1 relation between the rate constants in the model and the transformed parameters fitted in the model, the parameters obtained by the isometric logratio transformation are calculated from the set of formation fractions that quantify the paths to each of the compounds formed from a specific parent compound, and no such 1:1 relation exists.
+Therefore, parameter confidence intervals for formation fractions obtained with this method only appear valid for the case of a single transformation product, where only one formation fraction is to be estimated, directly corresponding to one component of the ilr transformed parameter.
+The confidence intervals obtained by backtransformation for the cases where a 1:1 relation between transformed and original parameter exist are considered by the author of this vignette to be more accurate than those obtained using a re-estimation of the Hessian matrix after backtransformation, as implemented in the FME package.
+
+
+Parameter t-test based on untransformed parameters
+The standard output of many nonlinear regression software packages includes the results from a test for significant difference from zero for all parameters. Such a test is also recommended to check the validity of rate constants in the FOCUS guidance (FOCUS Work Group on Degradation Kinetics 2014, 96ff).
+It has been argued that the precondition for this test, i.e. normal distribution of the estimator for the parameters, is not fulfilled in the case of nonlinear regression (Ranke and Lehmann 2015). However, this test is commonly used by industry, consultants and national authorities in order to decide on the reliability of parameter estimates, based on the FOCUS guidance mentioned above. Therefore, the results of this one-sided t-test are included in the summary output from mkin.
+As it is not reasonable to test for significant difference of the transformed parameters (e.g. \(log(k)\)) from zero, the t-test is calculated based on the model definition before parameter transformation, i.e. in a similar way as in packages that do not apply such an internal parameter transformation. A note is included in the mkin output, pointing to the fact that the t-test is based on the unjustified assumption of normal distribution of the parameter estimators.
+
+
+
+References
+
+
+Bates, D., and D. Watts. 1988. Nonlinear Regression and Its Applications. Wiley-Interscience.
+FOCUS Work Group on Degradation Kinetics. 2006. Guidance Document on Estimating Persistence and Degradation Kinetics from Environmental Fate Studies on Pesticides in EU Registration. Report of the FOCUS Work Group on Degradation Kinetics. http://focus.jrc.ec.europa.eu/dk.
+———. 2014. Generic Guidance for Estimating Persistence and Degradation Kinetics from Environmental Fate Studies on Pesticides in EU Registration. 1.1. http://focus.jrc.ec.europa.eu/dk.
+R Development Core Team. 2016. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org.
+Ranke, J. 2015. `kinfit`: Routines for Fitting Simple Kinetic Models to Chemical Degradation Data. http://CRAN.R-project.org/package=kinfit.
+———. 2016. `mkin`: Kinetic Evaluation of Chemical Degradation Data. http://CRAN.R-project.org/package=mkin.
+Ranke, J., and R. Lehmann. 2012. “Parameter Reliability in Kinetic Evaluation of Environmental Metabolism Data - Assessment and the Influence of Model Specification.” In SETAC World 20-24 May. Berlin.
+———. 2015. “To T-Test or Not to T-Test, That Is the Question.” In XV Symposium on Pesticide Chemistry 2-4 September 2015. Piacenza. http://chem.uft.uni-bremen.de/ranke/posters/piacenza_2015.pdf.
+Schäfer, D., B. Mikolasch, P. Rainbird, and B. Harvey. 2007. “KinGUI: a New Kinetic Software Tool for Evaluations According to FOCUS Degradation Kinetics.” In Proceedings of the XIII Symposium Pesticide Chemistry, edited by Del Re A. A. M., Capri E., Fragoulis G., and Trevisan M., 916–23. Piacenza.
+Soetaert, Karline, and Thomas Petzoldt. 2010. “Inverse Modelling, Sensitivity and Monte Carlo Analysis in R Using Package FME.” Journal of Statistical Software 33 (3): 1–28. http://www.jstatsoft.org/v33/i03/.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-- 
cgit v1.2.1

Example evaluation of FOCUS Example Dataset D

Example evaluation of FOCUS Example Dataset D

Johannes Ranke

2016-03-24

2016-06-28

Example evaluation of FOCUS Laboratory Data L1 to L3

Johannes Ranke

2016-03-24

Example evaluation of FOCUS Laboratory Data L1 to L3

Johannes Ranke

2016-06-28