From 1ef7008be2a72a0847064ad9c2ddcfa16b055482 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 3 May 2019 19:14:15 +0200 Subject: Improve error model fitting Now we have a three stage fitting process for nonconstant error models: - Unweighted least squares - Only optimize the error model - Optimize both Static documentation rebuilt by pkgdown --- R/mkinds.R | 1 - R/mkinfit.R | 50 ++- build.log | 4 +- check.log | 1 + docs/articles/FOCUS_D.html | 42 +-- docs/articles/FOCUS_D_files/figure-html/plot-1.png | Bin 98391 -> 98398 bytes docs/articles/FOCUS_L.html | 128 +++---- .../figure-html/unnamed-chunk-10-1.png | Bin 29150 -> 29148 bytes .../figure-html/unnamed-chunk-12-1.png | Bin 54997 -> 54995 bytes .../figure-html/unnamed-chunk-6-1.png | Bin 23974 -> 23962 bytes docs/articles/mkin.html | 2 +- docs/articles/twa.html | 2 +- docs/articles/web_only/FOCUS_Z.html | 138 ++++--- .../figure-html/FOCUS_2006_Z_fits_10-1.png | Bin 128677 -> 128674 bytes .../figure-html/FOCUS_2006_Z_fits_11-1.png | Bin 127690 -> 127674 bytes .../figure-html/FOCUS_2006_Z_fits_11a-1.png | Bin 96509 -> 96502 bytes .../figure-html/FOCUS_2006_Z_fits_11b-1.png | Bin 22308 -> 22316 bytes .../figure-html/FOCUS_2006_Z_fits_6-1.png | Bin 128965 -> 128979 bytes .../figure-html/FOCUS_2006_Z_fits_9-1.png | Bin 108421 -> 108437 bytes docs/articles/web_only/NAFTA_examples.html | 145 ++++---- .../NAFTA_examples_files/figure-html/p14-1.png | Bin 54038 -> 54036 bytes .../NAFTA_examples_files/figure-html/p9b-1.png | Bin 50012 -> 50009 bytes docs/articles/web_only/benchmarks.html | 112 +++--- docs/articles/web_only/compiled_models.html | 12 +- docs/reference/Extract.mmkin.html | 4 +- docs/reference/IORE.solution.html | 6 +- docs/reference/NAFTA_SOP_Attachment-1.png | Bin 42265 -> 42265 bytes docs/reference/NAFTA_SOP_Attachment.html | 4 +- docs/reference/add_err-1.png | Bin 89823 -> 89826 bytes docs/reference/logLik.mkinfit.html | 2 +- docs/reference/logistic.solution.html | 10 +- docs/reference/mccall81_245T.html | 42 +-- docs/reference/mkinfit.html | 402 +++++++++++---------- docs/reference/mkinmod.html | 2 +- docs/reference/mkinparplot-1.png | Bin 16555 -> 16549 bytes docs/reference/mkinparplot.html | 3 +- docs/reference/mkinpredict.html | 18 +- docs/reference/mkinresplot-1.png | Bin 14889 -> 14892 bytes docs/reference/mmkin-2.png | Bin 88227 -> 88226 bytes docs/reference/mmkin-3.png | Bin 85742 -> 85737 bytes docs/reference/mmkin.html | 10 +- docs/reference/plot.mkinfit-1.png | Bin 45173 -> 45182 bytes docs/reference/plot.mmkin-3.png | Bin 25441 -> 25445 bytes docs/reference/schaefer07_complex_case.html | 10 +- docs/reference/summary.mkinfit.html | 12 +- docs/reference/test_data_from_UBA_2014.html | 34 +- test.log | 16 +- tests/testthat/FOCUS_2006_D.csf | 2 +- vignettes/mkin_benchmarks.rda | Bin 801 -> 794 bytes 49 files changed, 622 insertions(+), 592 deletions(-) diff --git a/R/mkinds.R b/R/mkinds.R index 257a17c4..5333ca26 100644 --- a/R/mkinds.R +++ b/R/mkinds.R @@ -42,7 +42,6 @@ mkinds <- R6Class("mkinds", ) ) -#' @export print.mkinds <- function(x, ...) { cat(" with $title: ", x$title, "\n") cat("Observed compounds $observed: ", paste(x$observed, collapse = ", "), "\n") diff --git a/R/mkinfit.R b/R/mkinfit.R index dca71ecf..cca75690 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -263,7 +263,7 @@ mkinfit <- function(mkinmod, observed, errparms = rep(3, length(obs_vars)) } if (err_mod == "tc") { - errparms <- c(sigma_low = 3, rsd_high = 0.01) + errparms <- c(sigma_low = 0.1, rsd_high = 0.1) } names(errparms) <- errparm_names } @@ -275,13 +275,20 @@ mkinfit <- function(mkinmod, observed, length.out = n.outtimes)))) # Define log-likelihood function for optimisation, including (back)transformations - nlogLik <- function(P, trans = TRUE, OLS = FALSE, update_data = TRUE, ...) + nlogLik <- function(P, trans = TRUE, OLS = FALSE, local = FALSE, update_data = TRUE, ...) { assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter + P.orig <- P # Trace parameter values if requested if(trace_parms) cat(P, "\n") + # If we do a local optimisation of the error model, the initials + # for the state variabels and the parameters are given as 'local' + if (local[1] != FALSE) { + P <- local + } + # Initial states for t0 if(length(state.ini.optim) > 0) { odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed) @@ -291,7 +298,7 @@ mkinfit <- function(mkinmod, observed, names(odeini) <- state.ini.fixed.boxnames } - if (OLS) { + if (OLS | identical(P, local)) { odeparms.optim <- P[(length(state.ini.optim) + 1):length(P)] } else { odeparms.optim <- P[(length(state.ini.optim) + 1):(length(P) - length(errparms))] @@ -314,6 +321,9 @@ mkinfit <- function(mkinmod, observed, method.ode = method.ode, atol = atol, rtol = rtol, ...) + # Get back the original parameter vector to get the error model parameters + P <- P.orig + out_long <- mkin_wide_to_long(out, time = "time") if (err_mod == "const") { @@ -411,16 +421,17 @@ mkinfit <- function(mkinmod, observed, # Do the fit and take the time until the hessians are calculated fit_time <- system.time({ - # For constant variance, we do not include sigma in the optimisation, as it - # is unnecessary and leads to instability of the fit which are most obvious - # when fitting the IORE model - if (err_mod == "const") { - fit.ols <- nlminb(c(state.ini.optim, transparms.optim), - nlogLik, control = control, - lower = lower, upper = upper, OLS = TRUE, ...) + # In the inital run, we assume a constant standard deviation and do + # not optimise it, as this is unnecessary and leads to instability of the + # fit (most obvious when fitting the IORE model) + if (!quiet) message("Ordinary least squares optimisation") + parms.start <- c(state.ini.optim, transparms.optim) + fit.ols <- nlminb(parms.start, nlogLik, control = control, + lower = lower[names(parms.start)], + upper = upper[names(parms.start)], OLS = TRUE, ...) + if (err_mod == "const") { # Get the maximum likelihood estimate for sigma at the optimum parameter values - # This is equivalent to including sigma in the optimisation, but more stable data_errmod$residual <- data_errmod$value.observed - data_errmod$value.predicted sigma_mle = sqrt(sum(data_errmod$residual^2)/nrow(data_errmod)) @@ -428,9 +439,20 @@ mkinfit <- function(mkinmod, observed, fit <- fit.ols fit$logLik <- - nlogLik(c(fit$par, sigma = sigma_mle), OLS = FALSE) } else { - fit <- nlminb(c(state.ini.optim, transparms.optim, errparms), - nlogLik, control = control, - lower = lower, upper = upper, ...) + # After the OLS stop we fit the error model with fixed degradation model + # parameters + if (!quiet) message("Optimising the error model") + fit.err <- nlminb(errparms, nlogLik, control = control, + lower = lower[names(errparms)], + upper = upper[names(errparms)], + local = fit.ols$par, ...) + errparms.tmp <- fit.err$par + if (!quiet) message("Optimising the complete model") + parms.start <- c(state.ini.optim, transparms.optim, errparms.tmp) + fit <- nlminb(parms.start, nlogLik, + lower = lower[names(parms.start)], + upper = upper[names(parms.start)], + control = control, ...) fit$logLik <- - nlogLik.current } diff --git a/build.log b/build.log index 980d0de0..9728241a 100644 --- a/build.log +++ b/build.log @@ -6,7 +6,5 @@ * checking for LF line-endings in source and make files and shell scripts * checking for empty or unneeded directories * looking to see if a ‘data/datalist’ file should be added - NB: this package now depends on R (>= 3.5.0) - WARNING: Added dependency on R >= 3.5.0 because serialized objects in serialize/load version 3 cannot be read in older versions of R. File(s) containing such objects: 'mkin/vignettes/mkin_benchmarks.rda' -* building 'mkin_0.9.49.4.tar.gz' +* building ‘mkin_0.9.49.4.tar.gz’ diff --git a/check.log b/check.log index 325b690c..49ba0af6 100644 --- a/check.log +++ b/check.log @@ -17,6 +17,7 @@ Maintainer: ‘Johannes Ranke ’ * checking for hidden files and directories ... OK * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK +* checking serialization versions ... OK * checking whether package ‘mkin’ can be installed ... OK * checking installed package size ... OK * checking package directory ... OK diff --git a/docs/articles/FOCUS_D.html b/docs/articles/FOCUS_D.html index 3ce94de7..ae895cc3 100644 --- a/docs/articles/FOCUS_D.html +++ b/docs/articles/FOCUS_D.html @@ -88,7 +88,7 @@

Example evaluation of FOCUS Example Dataset D

Johannes Ranke

-

2019-05-02

+

2019-05-03

@@ -168,8 +168,8 @@ 
summary(fit)
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:13 2019 
-## Date of summary: Thu May  2 18:49:14 2019 
+## Date of fit:     Fri May  3 19:08:45 2019 
+## Date of summary: Fri May  3 19:08:46 2019 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent_sink * parent - k_parent_m1 * parent
@@ -177,7 +177,7 @@
 ## 
 ## Model predictions using solution type deSolve 
 ## 
-## Fitted using 396 model solutions performed in 0.999 s
+## Fitted using 389 model solutions performed in 0.983 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -216,13 +216,13 @@
 ## log_k_parent_sink  6.067e-01         1.000e+00      -8.550e-02
 ## log_k_parent_m1   -6.372e-02        -8.550e-02       1.000e+00
 ## log_k_m1_sink     -1.688e-01        -6.252e-01       4.731e-01
-## sigma             -3.368e-07         3.365e-08       8.420e-07
+## sigma              1.164e-09        -8.908e-10       1.652e-08
 ##                   log_k_m1_sink      sigma
-## parent_0             -1.688e-01 -3.368e-07
-## log_k_parent_sink    -6.252e-01  3.365e-08
-## log_k_parent_m1       4.731e-01  8.420e-07
-## log_k_m1_sink         1.000e+00  1.958e-08
-## sigma                 1.958e-08  1.000e+00
+## parent_0             -1.688e-01  1.164e-09
+## log_k_parent_sink    -6.252e-01 -8.908e-10
+## log_k_parent_m1       4.731e-01  1.652e-08
+## log_k_m1_sink         1.000e+00 -1.340e-10
+## sigma                -1.340e-10  1.000e+00
 ## 
 ## Backtransformed parameters:
 ## Confidence intervals for internally transformed parameters are asymmetric.
@@ -258,10 +258,10 @@
 ##     0   parent   102.04  99.59848  2.442e+00
 ##     1   parent    93.50  90.23787  3.262e+00
 ##     1   parent    92.50  90.23787  2.262e+00
-##     3   parent    63.23  74.07319 -1.084e+01
-##     3   parent    68.99  74.07319 -5.083e+00
-##     7   parent    52.32  49.91206  2.408e+00
-##     7   parent    55.13  49.91206  5.218e+00
+##     3   parent    63.23  74.07320 -1.084e+01
+##     3   parent    68.99  74.07320 -5.083e+00
+##     7   parent    52.32  49.91207  2.408e+00
+##     7   parent    55.13  49.91207  5.218e+00
 ##    14   parent    27.27  25.01257  2.257e+00
 ##    14   parent    26.64  25.01257  1.627e+00
 ##    21   parent    11.50  12.53462 -1.035e+00
@@ -271,23 +271,23 @@
 ##    50   parent     0.69   0.71624 -2.624e-02
 ##    50   parent     0.63   0.71624 -8.624e-02
 ##    75   parent     0.05   0.06074 -1.074e-02
-##    75   parent     0.06   0.06074 -7.381e-04
+##    75   parent     0.06   0.06074 -7.382e-04
 ##     1       m1     4.84   4.80296  3.704e-02
 ##     1       m1     5.64   4.80296  8.370e-01
 ##     3       m1    12.91  13.02400 -1.140e-01
 ##     3       m1    12.96  13.02400 -6.400e-02
 ##     7       m1    22.97  25.04476 -2.075e+00
 ##     7       m1    24.47  25.04476 -5.748e-01
-##    14       m1    41.69  36.69003  5.000e+00
-##    14       m1    33.21  36.69003 -3.480e+00
+##    14       m1    41.69  36.69002  5.000e+00
+##    14       m1    33.21  36.69002 -3.480e+00
 ##    21       m1    44.37  41.65310  2.717e+00
 ##    21       m1    46.44  41.65310  4.787e+00
 ##    35       m1    41.22  43.31312 -2.093e+00
 ##    35       m1    37.95  43.31312 -5.363e+00
-##    50       m1    41.19  41.21832 -2.832e-02
-##    50       m1    40.01  41.21832 -1.208e+00
-##    75       m1    40.09  36.44704  3.643e+00
-##    75       m1    33.85  36.44704 -2.597e+00
+##    50       m1    41.19  41.21831 -2.831e-02
+##    50       m1    40.01  41.21831 -1.208e+00
+##    75       m1    40.09  36.44703  3.643e+00
+##    75       m1    33.85  36.44703 -2.597e+00
 ##   100       m1    31.04  31.98163 -9.416e-01
 ##   100       m1    33.13  31.98163  1.148e+00
 ##   120       m1    25.15  28.78984 -3.640e+00
diff --git a/docs/articles/FOCUS_D_files/figure-html/plot-1.png b/docs/articles/FOCUS_D_files/figure-html/plot-1.png
index b8f2fe94..0cd6596d 100644
Binary files a/docs/articles/FOCUS_D_files/figure-html/plot-1.png and b/docs/articles/FOCUS_D_files/figure-html/plot-1.png differ
diff --git a/docs/articles/FOCUS_L.html b/docs/articles/FOCUS_L.html
index 42f2dc60..9507139c 100644
--- a/docs/articles/FOCUS_L.html
+++ b/docs/articles/FOCUS_L.html
@@ -88,7 +88,7 @@
       
Example evaluation of FOCUS Laboratory Data L1 to L3

                         
Johannes Ranke

             
-            
2019-05-02

+            
2019-05-03

       
       
       
@@ -114,15 +114,15 @@
 summary(m.L1.SFO)
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:16 2019 
-## Date of summary: Thu May  2 18:49:16 2019 
+## Date of fit:     Fri May  3 19:08:47 2019 
+## Date of summary: Fri May  3 19:08:47 2019 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 133 model solutions performed in 0.278 s
+## Fitted using 133 model solutions performed in 0.281 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -215,8 +215,8 @@
 ## finite result is doubtful
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:17 2019 
-## Date of summary: Thu May  2 18:49:17 2019 
+## Date of fit:     Fri May  3 19:08:50 2019 
+## Date of summary: Fri May  3 19:08:50 2019 
 ## 
 ## 
 ## Warning: Optimisation did not converge:
@@ -228,7 +228,7 @@
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 599 model solutions performed in 1.26 s
+## Fitted using 899 model solutions performed in 1.885 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -238,41 +238,41 @@
 ## parent_0 89.850000  state
 ## alpha     1.000000 deparm
 ## beta     10.000000 deparm
-## sigma     2.779868  error
+## sigma     2.779871  error
 ## 
 ## Starting values for the transformed parameters actually optimised:
 ##               value lower upper
 ## parent_0  89.850000  -Inf   Inf
 ## log_alpha  0.000000  -Inf   Inf
 ## log_beta   2.302585  -Inf   Inf
-## sigma      2.779868     0   Inf
+## sigma      2.779871     0   Inf
 ## 
 ## Fixed parameter values:
 ## None
 ## 
 ## Optimised, transformed parameters with symmetric confidence intervals:
 ##           Estimate Std. Error  Lower  Upper
-## parent_0     92.47     1.2810 89.720 95.220
-## log_alpha    10.66        NaN    NaN    NaN
-## log_beta     13.01        NaN    NaN    NaN
-## sigma         2.78     0.4599  1.794  3.766
+## parent_0     92.47     1.2800 89.730 95.220
+## log_alpha    10.58        NaN    NaN    NaN
+## log_beta     12.93        NaN    NaN    NaN
+## sigma         2.78     0.4507  1.813  3.747
 ## 
 ## Parameter correlation:
-##           parent_0 log_alpha log_beta    sigma
-## parent_0  1.000000       NaN      NaN 0.003475
-## log_alpha      NaN         1      NaN      NaN
-## log_beta       NaN       NaN        1      NaN
-## sigma     0.003475       NaN      NaN 1.000000
+##           parent_0 log_alpha log_beta   sigma
+## parent_0   1.00000       NaN      NaN 0.01452
+## log_alpha      NaN         1      NaN     NaN
+## log_beta       NaN       NaN        1     NaN
+## sigma      0.01452       NaN      NaN 1.00000
 ## 
 ## 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     92.47 72.13000 1.052e-19 89.720 95.220
-## alpha     42700.00  0.02298 4.910e-01     NA     NA
-## beta     446600.00  0.02298 4.910e-01     NA     NA
-## sigma         2.78  6.00000 1.628e-05  1.794  3.766
+## parent_0     92.47 72.13000 1.052e-19 89.730 95.220
+## alpha     39440.00  0.02397 4.906e-01     NA     NA
+## beta     412500.00  0.02397 4.906e-01     NA     NA
+## sigma         2.78  6.00000 1.628e-05  1.813  3.747
 ## 
 ## FOCUS Chi2 error levels in percent:
 ##          err.min n.optim df
@@ -281,7 +281,7 @@
 ## 
 ## Estimated disappearance times:
 ##         DT50  DT90 DT50back
-## parent 7.249 24.08     7.25
+## parent 7.249 24.08 7.249

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 \(\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 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 \(\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 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 \(\chi^2\) error levels was compared with KinGUII, CAKE and DegKin manager in a project sponsored by the German Umweltbundesamt (Ranke 2014).

@@ -319,15 +319,15 @@ 
summary(m.L2.FOMC, data = FALSE)
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:18 2019 
-## Date of summary: Thu May  2 18:49:18 2019 
+## Date of fit:     Fri May  3 19:08:51 2019 
+## Date of summary: Fri May  3 19:08:51 2019 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 240 model solutions performed in 0.486 s
+## Fitted using 239 model solutions performed in 0.486 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -358,10 +358,10 @@
 ## 
 ## Parameter correlation:
 ##             parent_0  log_alpha   log_beta      sigma
-## parent_0   1.000e+00 -1.151e-01 -2.085e-01  1.606e-08
-## log_alpha -1.151e-01  1.000e+00  9.741e-01 -1.168e-07
-## log_beta  -2.085e-01  9.741e-01  1.000e+00 -1.029e-07
-## sigma      1.606e-08 -1.168e-07 -1.029e-07  1.000e+00
+## parent_0   1.000e+00 -1.151e-01 -2.085e-01 -7.637e-09
+## log_alpha -1.151e-01  1.000e+00  9.741e-01 -1.617e-07
+## log_beta  -2.085e-01  9.741e-01  1.000e+00 -1.387e-07
+## sigma     -7.637e-09 -1.617e-07 -1.387e-07  1.000e+00
 ## 
 ## Backtransformed parameters:
 ## Confidence intervals for internally transformed parameters are asymmetric.
@@ -394,8 +394,8 @@
 
summary(m.L2.DFOP, data = FALSE)

 
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:20 2019 
-## Date of summary: Thu May  2 18:49:20 2019 
+## Date of fit:     Fri May  3 19:08:52 2019 
+## Date of summary: Fri May  3 19:08:52 2019 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) *
@@ -404,7 +404,7 @@
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 587 model solutions performed in 1.223 s
+## Fitted using 572 model solutions performed in 1.193 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -431,18 +431,18 @@
 ## Optimised, transformed parameters with symmetric confidence intervals:
 ##          Estimate Std. Error      Lower     Upper
 ## parent_0  93.9500  9.998e-01    91.5900   96.3100
-## log_k1     3.1330  2.265e+03 -5354.0000 5360.0000
+## log_k1     3.1370  2.376e+03 -5616.0000 5622.0000
 ## log_k2    -1.0880  6.285e-02    -1.2370   -0.9394
 ## g_ilr     -0.2821  7.033e-02    -0.4484   -0.1158
 ## sigma      1.4140  2.886e-01     0.7314    2.0960
 ## 
 ## Parameter correlation:
 ##            parent_0     log_k1     log_k2      g_ilr      sigma
-## parent_0  1.000e+00  5.434e-07 -9.989e-11  2.665e-01 -3.978e-10
-## log_k1    5.434e-07  1.000e+00  8.888e-05 -1.748e-04 -8.207e-06
-## log_k2   -9.989e-11  8.888e-05  1.000e+00 -7.903e-01  5.751e-10
-## g_ilr     2.665e-01 -1.748e-04 -7.903e-01  1.000e+00 -7.109e-10
-## sigma    -3.978e-10 -8.207e-06  5.751e-10 -7.109e-10  1.000e+00
+## parent_0  1.000e+00  5.155e-07  2.371e-09  2.665e-01 -6.849e-09
+## log_k1    5.155e-07  1.000e+00  8.434e-05 -1.659e-04 -7.791e-06
+## log_k2    2.371e-09  8.434e-05  1.000e+00 -7.903e-01 -1.262e-08
+## g_ilr     2.665e-01 -1.659e-04 -7.903e-01  1.000e+00  3.241e-08
+## sigma    -6.849e-09 -7.791e-06 -1.262e-08  3.241e-08  1.000e+00
 ## 
 ## Backtransformed parameters:
 ## Confidence intervals for internally transformed parameters are asymmetric.
@@ -450,7 +450,7 @@
 ## for estimators of untransformed parameters.
 ##          Estimate   t value    Pr(>t)   Lower   Upper
 ## parent_0  93.9500 9.397e+01 2.036e-12 91.5900 96.3100
-## k1        22.9300 4.514e-04 4.998e-01  0.0000     Inf
+## k1        23.0400 4.303e-04 4.998e-01  0.0000     Inf
 ## k2         0.3369 1.591e+01 4.697e-07  0.2904  0.3909
 ## g          0.4016 1.680e+01 3.238e-07  0.3466  0.4591
 ## sigma      1.4140 4.899e+00 8.776e-04  0.7314  2.0960
@@ -462,7 +462,7 @@
 ## 
 ## Estimated disappearance times:
 ##          DT50  DT90 DT50_k1 DT50_k2
-## parent 0.5335 5.311 0.03023   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.

@@ -493,8 +493,8 @@ 
summary(mm.L3[["DFOP", 1]])
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:21 2019 
-## Date of summary: Thu May  2 18:49:22 2019 
+## Date of fit:     Fri May  3 19:08:54 2019 
+## Date of summary: Fri May  3 19:08:54 2019 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) *
@@ -503,7 +503,7 @@
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 372 model solutions performed in 0.768 s
+## Fitted using 373 model solutions performed in 0.775 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -536,12 +536,12 @@
 ## sigma      1.0170    0.25430  0.2079   1.827000
 ## 
 ## Parameter correlation:
-##           parent_0     log_k1     log_k2      g_ilr      sigma
-## parent_0 1.000e+00  1.732e-01  2.282e-02  4.009e-01  1.660e-07
-## log_k1   1.732e-01  1.000e+00  4.945e-01 -5.809e-01  6.635e-08
-## log_k2   2.282e-02  4.945e-01  1.000e+00 -6.812e-01  3.880e-07
-## g_ilr    4.009e-01 -5.809e-01 -6.812e-01  1.000e+00 -3.822e-07
-## sigma    1.660e-07  6.635e-08  3.880e-07 -3.822e-07  1.000e+00
+##            parent_0     log_k1     log_k2      g_ilr      sigma
+## parent_0  1.000e+00  1.732e-01  2.282e-02  4.009e-01 -6.872e-07
+## log_k1    1.732e-01  1.000e+00  4.945e-01 -5.809e-01  3.200e-07
+## log_k2    2.282e-02  4.945e-01  1.000e+00 -6.812e-01  7.673e-07
+## g_ilr     4.009e-01 -5.809e-01 -6.812e-01  1.000e+00 -8.731e-07
+## sigma    -6.872e-07  3.200e-07  7.673e-07 -8.731e-07  1.000e+00
 ## 
 ## Backtransformed parameters:
 ## Confidence intervals for internally transformed parameters are asymmetric.
@@ -598,15 +598,15 @@
 
summary(mm.L4[["SFO", 1]], data = FALSE)

 
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:22 2019 
-## Date of summary: Thu May  2 18:49:23 2019 
+## Date of fit:     Fri May  3 19:08:55 2019 
+## Date of summary: Fri May  3 19:08:55 2019 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent_sink * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 146 model solutions performed in 0.294 s
+## Fitted using 142 model solutions performed in 0.29 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -633,10 +633,10 @@
 ## sigma                3.162    0.79050  1.130   5.194
 ## 
 ## Parameter correlation:
-##                    parent_0 log_k_parent_sink      sigma
-## parent_0          1.000e+00         5.938e-01  4.256e-10
-## log_k_parent_sink 5.938e-01         1.000e+00 -7.280e-10
-## sigma             4.256e-10        -7.280e-10  1.000e+00
+##                    parent_0 log_k_parent_sink     sigma
+## parent_0          1.000e+00         5.938e-01 3.440e-07
+## log_k_parent_sink 5.938e-01         1.000e+00 5.885e-07
+## sigma             3.440e-07         5.885e-07 1.000e+00
 ## 
 ## Backtransformed parameters:
 ## Confidence intervals for internally transformed parameters are asymmetric.
@@ -662,15 +662,15 @@
 
summary(mm.L4[["FOMC", 1]], data = FALSE)

 
## mkin version used for fitting:    0.9.49.4 
 ## R version used for fitting:       3.6.0 
-## Date of fit:     Thu May  2 18:49:23 2019 
-## Date of summary: Thu May  2 18:49:23 2019 
+## Date of fit:     Fri May  3 19:08:55 2019 
+## Date of summary: Fri May  3 19:08:55 2019 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 224 model solutions performed in 0.454 s
+## Fitted using 224 model solutions performed in 0.451 s
 ## 
 ## Error model:
 ## Constant variance 
@@ -701,10 +701,10 @@
 ## 
 ## Parameter correlation:
 ##             parent_0  log_alpha   log_beta      sigma
-## parent_0   1.000e+00 -4.696e-01 -5.543e-01 -2.473e-07
-## log_alpha -4.696e-01  1.000e+00  9.889e-01  2.429e-08
-## log_beta  -5.543e-01  9.889e-01  1.000e+00  5.183e-08
-## sigma     -2.473e-07  2.429e-08  5.183e-08  1.000e+00
+## parent_0   1.000e+00 -4.696e-01 -5.543e-01 -2.563e-07
+## log_alpha -4.696e-01  1.000e+00  9.889e-01  4.066e-08
+## log_beta  -5.543e-01  9.889e-01  1.000e+00  6.818e-08
+## sigma     -2.563e-07  4.066e-08  6.818e-08  1.000e+00
 ## 
 ## Backtransformed parameters:
 ## Confidence intervals for internally transformed parameters are asymmetric.
diff --git a/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-10-1.png
index 7f695ecf..cc28c092 100644
Binary files a/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-10-1.png and b/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-12-1.png b/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-12-1.png
index 850df1e1..46d85b73 100644
Binary files a/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-12-1.png and b/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-12-1.png differ
diff --git a/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-6-1.png b/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-6-1.png
index 5ed0b4de..3dd8565c 100644
Binary files a/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-6-1.png and b/docs/articles/FOCUS_L_files/figure-html/unnamed-chunk-6-1.png differ
diff --git a/docs/articles/mkin.html b/docs/articles/mkin.html
index 132e0c14..b0af83ce 100644
--- a/docs/articles/mkin.html
+++ b/docs/articles/mkin.html
@@ -88,7 +88,7 @@
       
Introduction to mkin

                         
Johannes Ranke

             
-            
2019-05-02

+            
2019-05-03

       
       
       
diff --git a/docs/articles/twa.html b/docs/articles/twa.html
index 7d9d1c4e..7f68d39e 100644
--- a/docs/articles/twa.html
+++ b/docs/articles/twa.html
@@ -88,7 +88,7 @@
       
Calculation of time weighted average concentrations with mkin

                         
Johannes Ranke

             
-            
2019-05-02

+            
2019-05-03

       
       
       
diff --git a/docs/articles/web_only/FOCUS_Z.html b/docs/articles/web_only/FOCUS_Z.html
index 950e8eab..9e64ae3a 100644
--- a/docs/articles/web_only/FOCUS_Z.html
+++ b/docs/articles/web_only/FOCUS_Z.html
@@ -88,7 +88,7 @@
       
Example evaluation of FOCUS dataset Z

                         
Johannes Ranke

             
-            
2019-05-02

+            
2019-05-03

       
       
       
@@ -132,11 +132,11 @@
 

 
summary(m.Z.2a, data = FALSE)$bpar

 
##             Estimate se_notrans    t value     Pr(>t)    Lower    Upper
-## Z0_0      9.7015e+01   3.394776 2.8578e+01 6.5093e-21 91.66556 102.3642
-## k_Z0_sink 4.0301e-10   0.225510 1.7871e-09 5.0000e-01  0.00000      Inf
-## k_Z0_Z1   2.2360e+00   0.159161 1.4049e+01 1.1412e-13  1.95303   2.5600
-## k_Z1_sink 4.8212e-01   0.065499 7.3608e+00 5.1791e-08  0.40341   0.5762
-## sigma     4.8041e+00   0.637657 7.5340e+00 3.4468e-08  3.52677   6.0815

+## Z0_0      9.7015e+01   3.393176 2.8591e+01 6.4352e-21 91.66556 102.3642
+## k_Z0_sink 7.2231e-10   0.225254 3.2067e-09 5.0000e-01  0.00000      Inf
+## k_Z0_Z1   2.2360e+00   0.159134 1.4051e+01 1.1369e-13  1.95303   2.5600
+## k_Z1_sink 4.8212e-01   0.065454 7.3658e+00 5.1186e-08  0.40341   0.5762
+## sigma     4.8041e+00   0.637618 7.5345e+00 3.4431e-08  3.52677   6.0815

 
As obvious from the parameter summary (the  component of the summary), the kinetic rate constant from parent compound Z to sink is very small and the t-test for this parameter suggests that it is not significantly different from zero. This suggests, in agreement with the analysis in the FOCUS kinetics report, to simplify the model by removing the pathway to sink.

 
A similar result can be obtained when formation fractions are used in the model formulation:

 
Z.2a.ff <- mkinmod(Z0 = mkinsub("SFO", "Z1"),
@@ -199,20 +199,18 @@
                      quiet = TRUE)

 
## Warning in mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, parms.ini = m.Z.
 ## 5$bparms.ode, : Observations with value of zero were removed from the data

-
## Warning in mkinfit(Z.FOCUS, FOCUS_2006_Z_mkin, parms.ini = m.Z.5$bparms.ode, : Optimisation did not converge:
-## false convergence (8)

-
plot_sep(m.Z.FOCUS)

+
plot_sep(m.Z.FOCUS)

 

-
summary(m.Z.FOCUS, data = FALSE)$bpar

+
summary(m.Z.FOCUS, data = FALSE)$bpar

 
##             Estimate se_notrans t value     Pr(>t)     Lower      Upper
-## Z0_0       96.838619   1.994272 48.5584 4.0282e-42 92.826596 100.850642
-## k_Z0        2.215408   0.118459 18.7018 1.0415e-23  1.989468   2.467007
-## k_Z1        0.478300   0.028257 16.9267 6.2407e-22  0.424701   0.538663
-## k_Z2        0.451618   0.042138 10.7177 1.6308e-14  0.374327   0.544869
-## k_Z3        0.058693   0.015246  3.8498 1.7805e-04  0.034804   0.098981
-## f_Z2_to_Z3  0.471508   0.058352  8.0804 9.6647e-11  0.357725   0.588332
+## Z0_0       96.838607   1.994273 48.5584 4.0283e-42 92.826626 100.850589
+## k_Z0        2.215405   0.118459 18.7018 1.0415e-23  1.989465   2.467003
+## k_Z1        0.478300   0.028257 16.9267 6.2408e-22  0.424701   0.538662
+## k_Z2        0.451618   0.042138 10.7177 1.6308e-14  0.374328   0.544867
+## k_Z3        0.058693   0.015246  3.8498 1.7806e-04  0.034805   0.098978
+## f_Z2_to_Z3  0.471508   0.058352  8.0804 9.6648e-11  0.357735   0.588320
 ## sigma       3.984431   0.383402 10.3923 4.5575e-14  3.213126   4.755736

-
endpoints(m.Z.FOCUS)

+
endpoints(m.Z.FOCUS)

 
## $ff
 ##   Z2_Z3 Z2_sink 
 ## 0.47151 0.52849 
@@ -225,7 +223,7 @@
 ## Z0  0.31288  1.0394
 ## Z1  1.44919  4.8141
 ## Z2  1.53481  5.0985
-## Z3 11.80962 39.2307

+## Z3 11.80965 39.2308

 
This fit corresponds to the final result chosen in Appendix 7 of the FOCUS report. Confidence intervals returned by mkin are based on internally transformed parameters, however.

 
 

@@ -233,17 +231,17 @@
 Using the SFORB model
 
As the FOCUS report states, there is a certain tailing of the time course of metabolite Z3. Also, the time course of the parent compound is not fitted very well using the SFO model, as residues at a certain low level remain.

 
Therefore, an additional model is offered here, using the single first-order reversible binding (SFORB) model for metabolite Z3. As expected, the \(\chi^2\) error level is lower for metabolite Z3 using this model and the graphical fit for Z3 is improved. However, the covariance matrix is not returned.

-
Z.mkin.1 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE),
-                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
-                    Z2 = mkinsub("SFO", "Z3"),
-                    Z3 = mkinsub("SFORB"))

+
Z.mkin.1 <- mkinmod(Z0 = mkinsub("SFO", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO", "Z3"),
+                    Z3 = mkinsub("SFORB"))

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

-
m.Z.mkin.1 <- mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, quiet = TRUE)

+
m.Z.mkin.1 <- mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, quiet = TRUE)

 
## Warning in mkinfit(Z.mkin.1, FOCUS_2006_Z_mkin, quiet = TRUE): Observations
 ## with value of zero were removed from the data

-
plot_sep(m.Z.mkin.1)

+
plot_sep(m.Z.mkin.1)

 

-
summary(m.Z.mkin.1, data = FALSE)$cov.unscaled

+
summary(m.Z.mkin.1, data = FALSE)$cov.unscaled

 
##                            Z0_0 log_k_Z0_Z1 log_k_Z1_Z2 log_k_Z2_sink
 ## Z0_0                 3.8375e+00  5.4918e-03  3.0584e-02    1.2969e-01
 ## log_k_Z0_Z1          5.4918e-03  2.7613e-03 -1.8820e-04    2.6634e-04
@@ -251,9 +249,9 @@
 ## log_k_Z2_sink        1.2969e-01  2.6634e-04  3.2177e-03    3.4256e-02
 ## log_k_Z2_Z3_free    -2.4223e-02 -2.6169e-04 -1.1845e-03   -8.1134e-03
 ## log_k_Z3_free_sink  -6.5467e-02 -4.0815e-04 -3.2978e-03   -3.6010e-02
-## log_k_Z3_free_bound -6.0658e-02 -4.4768e-04 -3.0588e-03   -3.9074e-02
-## log_k_Z3_bound_free  4.7821e+00  5.5819e-03  1.0267e-01    1.1956e+00
-## sigma               -1.4345e-08  8.6519e-11 -6.1861e-10   -4.7499e-10
+## log_k_Z3_free_bound -6.0659e-02 -4.4768e-04 -3.0588e-03   -3.9074e-02
+## log_k_Z3_bound_free  5.2844e-01  4.5458e-03  7.9800e-03    4.6274e-02
+## sigma                2.0366e-10 -3.4658e-10  8.9910e-11   -2.5946e-10
 ##                     log_k_Z2_Z3_free log_k_Z3_free_sink
 ## Z0_0                     -2.4223e-02        -6.5467e-02
 ## log_k_Z0_Z1              -2.6169e-04        -4.0815e-04
@@ -262,84 +260,84 @@
 ## log_k_Z2_Z3_free          1.5500e-02         2.1583e-02
 ## log_k_Z3_free_sink        2.1583e-02         7.5705e-02
 ## log_k_Z3_free_bound       2.5836e-02         1.1964e-01
-## log_k_Z3_bound_free      -2.1303e-01        -9.0584e-01
-## sigma                     5.8776e-10         1.0773e-09
+## log_k_Z3_bound_free       5.2534e-02         2.9441e-01
+## sigma                     1.3063e-10         3.4170e-10
 ##                     log_k_Z3_free_bound log_k_Z3_bound_free       sigma
-## Z0_0                        -6.0658e-02          4.7821e+00 -1.4345e-08
-## log_k_Z0_Z1                 -4.4768e-04          5.5819e-03  8.6519e-11
-## log_k_Z1_Z2                 -3.0588e-03          1.0267e-01 -6.1861e-10
-## log_k_Z2_sink               -3.9074e-02          1.1956e+00 -4.7499e-10
-## log_k_Z2_Z3_free             2.5836e-02         -2.1303e-01  5.8776e-10
-## log_k_Z3_free_sink           1.1964e-01         -9.0584e-01  1.0773e-09
-## log_k_Z3_free_bound          6.5902e-01          4.2011e+00  2.1743e-09
-## log_k_Z3_bound_free          4.2011e+00          3.6036e+08  7.2404e-02
-## sigma                        2.1743e-09          7.2404e-02  1.4170e-01

+## Z0_0                        -6.0659e-02          5.2844e-01  2.0366e-10
+## log_k_Z0_Z1                 -4.4768e-04          4.5458e-03 -3.4658e-10
+## log_k_Z1_Z2                 -3.0588e-03          7.9800e-03  8.9910e-11
+## log_k_Z2_sink               -3.9074e-02          4.6274e-02 -2.5946e-10
+## log_k_Z2_Z3_free             2.5836e-02          5.2534e-02  1.3063e-10
+## log_k_Z3_free_sink           1.1964e-01          2.9441e-01  3.4170e-10
+## log_k_Z3_free_bound          6.5902e-01          5.4737e+00 -6.7704e-10
+## log_k_Z3_bound_free          5.4737e+00          2.8722e+08  7.2421e-02
+## sigma                       -6.7704e-10          7.2421e-02  1.4170e-01

Therefore, a further stepwise model building is performed starting from the stage of parent and two metabolites, starting from the assumption that the model fit for the parent compound can be improved by using the SFORB model.

-
Z.mkin.3 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
-                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
-                    Z2 = mkinsub("SFO"))
+
Z.mkin.3 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO"))
## Successfully compiled differential equation model from auto-generated C code.
-
m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE)
+
m.Z.mkin.3 <- mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE)
## Warning in mkinfit(Z.mkin.3, FOCUS_2006_Z_mkin, quiet = TRUE): Observations
 ## with value of zero were removed from the data
-
plot_sep(m.Z.mkin.3)
+
plot_sep(m.Z.mkin.3)

This results in a much better representation of the behaviour of the parent compound Z0.

Finally, Z3 is added as well. These models appear overparameterised (no covariance matrix returned) if the sink for Z1 is left in the models.

-
Z.mkin.4 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
-                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
-                    Z2 = mkinsub("SFO", "Z3"),
-                    Z3 = mkinsub("SFO"))
+
Z.mkin.4 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO", "Z3"),
+                    Z3 = mkinsub("SFO"))
## Successfully compiled differential equation model from auto-generated C code.
- + 
## Warning in mkinfit(Z.mkin.4, FOCUS_2006_Z_mkin, parms.ini = m.Z.mkin.
 ## 3$bparms.ode, : Observations with value of zero were removed from the data
-
plot_sep(m.Z.mkin.4)
+
plot_sep(m.Z.mkin.4)

The error level of the fit, but especially of metabolite Z3, can be improved if the SFORB model is chosen for this metabolite, as this model is capable of representing the tailing of the metabolite decline phase.

-
Z.mkin.5 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
-                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
-                    Z2 = mkinsub("SFO", "Z3"),
-                    Z3 = mkinsub("SFORB"))
+
Z.mkin.5 <- mkinmod(Z0 = mkinsub("SFORB", "Z1", sink = FALSE),
+                    Z1 = mkinsub("SFO", "Z2", sink = FALSE),
+                    Z2 = mkinsub("SFO", "Z3"),
+                    Z3 = mkinsub("SFORB"))
## Successfully compiled differential equation model from auto-generated C code.
- + 
## Warning in mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, parms.ini = m.Z.mkin.
 ## 4$bparms.ode[1:4], : Observations with value of zero were removed from the
 ## data
-
plot_sep(m.Z.mkin.5)
+
plot_sep(m.Z.mkin.5)

The summary view of the backtransformed parameters shows that we get no confidence intervals due to overparameterisation. As the optimized is excessively small, it seems reasonable to fix it to zero.

- + 
## Warning in mkinfit(Z.mkin.5, FOCUS_2006_Z_mkin, parms.ini = c(m.Z.mkin.
 ## 5$bparms.ode[1:7], : Observations with value of zero were removed from the
 ## data
-
plot_sep(m.Z.mkin.5a)
+
plot_sep(m.Z.mkin.5a)

As expected, the residual plots for Z0 and Z3 are more random than in the case of the all SFO model for which they were shown above. In conclusion, the model is proposed as the best-fit model for the dataset from Appendix 7 of the FOCUS report.

A graphical representation of the confidence intervals can finally be obtained.

-
mkinparplot(m.Z.mkin.5a)
+
mkinparplot(m.Z.mkin.5a)

The endpoints obtained with this model are

-
endpoints(m.Z.mkin.5a)
+
endpoints(m.Z.mkin.5a)
## $ff
 ##   Z0_free_Z1        Z1_Z2      Z2_sink   Z2_Z3_free Z3_free_sink 
 ##      1.00000      1.00000      0.46344      0.53656      1.00000 
 ## 
 ## $SFORB
 ##     Z0_b1     Z0_b2     Z3_b1     Z3_b2 
-## 2.4471329 0.0075123 0.0800074 0.0000000 
+## 2.4471381 0.0075124 0.0800075 0.0000000 
 ## 
 ## $distimes
 ##      DT50   DT90 DT50_Z0_b1 DT50_Z0_b2 DT50_Z3_b1 DT50_Z3_b2
-## Z0 0.3043 1.1848    0.28325     92.268         NA         NA
+## Z0 0.3043 1.1848    0.28325     92.267         NA         NA
 ## Z1 1.5148 5.0320         NA         NA         NA         NA
 ## Z2 1.6414 5.4526         NA         NA         NA         NA
 ## Z3     NA     NA         NA         NA     8.6635        Inf
diff --git a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png index 3dfb6477..1a92365e 100644 Binary files a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png and b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_10-1.png differ diff --git a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png index bb909e98..debf3cb8 100644 Binary files a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png and b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11-1.png differ diff --git a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png index c7cbb448..89020a0a 100644 Binary files a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png and b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11a-1.png differ diff --git a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png index e39a2e6a..8188c4cf 100644 Binary files a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png and b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_11b-1.png differ diff --git a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png index 7213dac1..a11e0092 100644 Binary files a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png and b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_6-1.png differ diff --git a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png index 27c51525..acd197f2 100644 Binary files a/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png and b/docs/articles/web_only/FOCUS_Z_files/figure-html/FOCUS_2006_Z_fits_9-1.png differ diff --git a/docs/articles/web_only/NAFTA_examples.html b/docs/articles/web_only/NAFTA_examples.html index 0ae02f3c..04a60db6 100644 --- a/docs/articles/web_only/NAFTA_examples.html +++ b/docs/articles/web_only/NAFTA_examples.html @@ -88,7 +88,7 @@

Evaluation of example datasets from Attachment 1 to the US EPA SOP for the NAFTA guidance

Johannes Ranke

-

2019-05-02

+

2019-05-03

@@ -141,7 +141,7 @@ ## Estimate Pr(>t) Lower Upper ## parent_0 9.99e+01 1.41e-26 98.8116 101.0810 ## k1 2.67e-02 5.05e-06 0.0243 0.0295 -## k2 3.41e-12 5.00e-01 0.0000 Inf +## k2 2.86e-12 5.00e-01 0.0000 Inf ## g 6.47e-01 3.67e-06 0.6248 0.6677 ## sigma 1.27e+00 8.91e-06 0.8395 1.6929 ## @@ -150,7 +150,7 @@ ## DT50 DT90 DT50_rep ## SFO 67.7 2.25e+02 6.77e+01 ## IORE 58.2 1.07e+03 3.22e+02 -## DFOP 55.5 3.70e+11 2.03e+11 +## DFOP 55.5 4.42e+11 2.42e+11 ## ## Representative half-life: ## [1] 321.51 @@ -189,7 +189,7 @@ ## Estimate Pr(>t) Lower Upper ## parent_0 9.84e+01 1.24e-27 97.8078 98.9187 ## k1 1.55e-02 4.10e-04 0.0143 0.0167 -## k2 8.28e-12 5.00e-01 0.0000 Inf +## k2 1.16e-11 5.00e-01 0.0000 Inf ## g 6.89e-01 2.92e-03 0.6626 0.7142 ## sigma 6.48e-01 2.38e-05 0.4147 0.8813 ## @@ -198,7 +198,7 @@ ## DT50 DT90 DT50_rep ## SFO 86.6 2.88e+02 8.66e+01 ## IORE 85.5 7.17e+02 2.16e+02 -## DFOP 83.6 1.37e+11 8.37e+10 +## DFOP 83.6 9.80e+10 5.98e+10 ## ## Representative half-life: ## [1] 215.87 @@ -222,7 +222,7 @@ ## Parameters: ## $SFO ## Estimate Pr(>t) Lower Upper -## parent_0 94.7759 7.29e-24 92.3479 97.2039 +## parent_0 94.7759 7.29e-24 92.3478 97.2039 ## k_parent_sink 0.0179 8.02e-16 0.0166 0.0194 ## sigma 3.0696 3.81e-06 2.0456 4.0936 ## @@ -237,7 +237,7 @@ ## Estimate Pr(>t) Lower Upper ## parent_0 9.66e+01 1.57e-25 95.3476 97.8979 ## k1 2.55e-02 7.33e-06 0.0233 0.0278 -## k2 3.87e-11 5.00e-01 0.0000 Inf +## k2 4.90e-11 5.00e-01 0.0000 Inf ## g 8.61e-01 7.55e-06 0.8314 0.8867 ## sigma 1.46e+00 6.93e-06 0.9661 1.9483 ## @@ -246,7 +246,7 @@ ## DT50 DT90 DT50_rep ## SFO 38.6 1.28e+02 3.86e+01 ## IORE 34.0 1.77e+02 5.32e+01 -## DFOP 34.1 8.43e+09 1.79e+10 +## DFOP 34.1 6.66e+09 1.41e+10 ## ## Representative half-life: ## [1] 53.17 @@ -285,7 +285,7 @@ ## Estimate Pr(>t) Lower Upper ## parent_0 9.89e+01 9.44e-49 95.4640 102.2573 ## k1 1.81e-02 1.75e-01 0.0116 0.0281 -## k2 2.46e-10 5.00e-01 0.0000 Inf +## k2 1.97e-10 5.00e-01 0.0000 Inf ## g 6.06e-01 2.19e-01 0.4826 0.7178 ## sigma 7.40e+00 2.97e-15 6.0201 8.7754 ## @@ -294,7 +294,7 @@ ## DT50 DT90 DT50_rep ## SFO 94.3 3.13e+02 9.43e+01 ## IORE 96.7 1.51e+03 4.55e+02 -## DFOP 96.4 5.58e+09 2.82e+09 +## DFOP 96.4 6.97e+09 3.52e+09 ## ## Representative half-life: ## [1] 454.55 @@ -397,7 +397,7 @@ ## Estimate Pr(>t) Lower Upper ## parent_0 9.85e+01 2.54e-20 97.390 99.672 ## k1 1.38e-01 3.52e-05 0.131 0.146 -## k2 5.75e-13 5.00e-01 0.000 Inf +## k2 6.02e-13 5.00e-01 0.000 Inf ## g 6.52e-01 8.13e-06 0.642 0.661 ## sigma 7.88e-01 6.13e-02 0.481 1.095 ## @@ -406,7 +406,7 @@ ## DT50 DT90 DT50_rep ## SFO 16.9 5.63e+01 1.69e+01 ## IORE 11.6 3.37e+02 1.01e+02 -## DFOP 10.5 2.17e+12 1.21e+12 +## DFOP 10.5 2.07e+12 1.15e+12 ## ## Representative half-life: ## [1] 101.43 @@ -416,12 +416,16 @@

Example on page 9, lower panel

p9b <- nafta(NAFTA_SOP_Attachment[["p9b"]])
+
## Warning in sqrt(diag(covar)): NaNs wurden erzeugt
## Warning in sqrt(diag(covar_notrans)): NaNs wurden erzeugt
+
## Warning in sqrt(1/diag(V)): NaNs wurden erzeugt
+
## Warning in cov2cor(ans$cov.unscaled): diag(.) had 0 or NA entries; non-
+## finite result is doubtful
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c
## The half-life obtained from the IORE model may be used
-
plot(p9b)
+
plot(p9b)

-
print(p9b)
+
print(p9b)
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 35.64867 23.22334 35.64867 
@@ -446,9 +450,9 @@
 ## $DFOP
 ##          Estimate   Pr(>t)   Lower   Upper
 ## parent_0  94.7123 1.61e-16 93.1355 96.2891
-## k1         0.0389      NaN  0.0306  0.0495
-## k2         0.0389 1.10e-06  0.0186  0.0812
-## g          0.7598      NaN  0.0000  1.0000
+## k1         0.0389 1.43e-06  0.0312  0.0485
+## k2         0.0389 6.67e-03  0.0186  0.0812
+## g          0.7742      NaN      NA      NA
 ## sigma      1.5957 2.50e-04  0.9135  2.2779
 ## 
 ## 
@@ -465,12 +469,12 @@
 

 

 Example on page 10

-
p10 <- nafta(NAFTA_SOP_Attachment[["p10"]])

+
p10 <- nafta(NAFTA_SOP_Attachment[["p10"]])

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p10)

+
plot(p10)

 

-
print(p10)

+
print(p10)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 899.4089 336.4348 899.4089 
@@ -495,8 +499,8 @@
 ## $DFOP
 ##          Estimate   Pr(>t)   Lower    Upper
 ## parent_0 101.7315 1.41e-09 91.6534 111.8097
-## k1         0.0495 6.41e-04  0.0303   0.0809
-## k2         0.0495 1.66e-02  0.0201   0.1219
+## k1         0.0495 6.42e-04  0.0301   0.0814
+## k2         0.0495 1.66e-02  0.0200   0.1225
 ## g          0.6634 5.00e-01  0.0000   1.0000
 ## sigma      8.0152 2.50e-04  4.5886  11.4418
 ## 
@@ -518,12 +522,12 @@
 

 

 Example on page 11

-
p11 <- nafta(NAFTA_SOP_Attachment[["p11"]])

+
p11 <- nafta(NAFTA_SOP_Attachment[["p11"]])

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p11)

+
plot(p11)

 

-
print(p11)

+
print(p11)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 579.6805 204.7932 144.7783 
@@ -549,7 +553,7 @@
 ##          Estimate   Pr(>t)    Lower    Upper
 ## parent_0 1.05e+02 9.47e-13  99.9990 109.1224
 ## k1       4.41e-02 5.95e-03   0.0296   0.0658
-## k2       9.20e-13 5.00e-01   0.0000      Inf
+## k2       7.25e-13 5.00e-01   0.0000      Inf
 ## g        3.22e-01 1.45e-03   0.2814   0.3650
 ## sigma    3.22e+00 3.52e-04   1.8410   4.5906
 ## 
@@ -558,7 +562,7 @@
 ##          DT50     DT90 DT50_rep
 ## SFO  2.16e+02 7.18e+02 2.16e+02
 ## IORE 9.73e+02 1.37e+08 4.11e+07
-## DFOP 3.31e+11 2.08e+12 7.53e+11
+## DFOP 4.21e+11 2.64e+12 9.56e+11
 ## 
 ## Representative half-life:
 ## [1] 41148169

@@ -572,19 +576,14 @@
 

 

 Example on page 12, upper panel

-
p12a <- nafta(NAFTA_SOP_Attachment[["p12a"]])

+
p12a <- nafta(NAFTA_SOP_Attachment[["p12a"]])

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

-
## Warning in sqrt(diag(covar)): NaNs wurden erzeugt

-
## Warning in sqrt(diag(covar_notrans)): NaNs wurden erzeugt

-
## Warning in sqrt(1/diag(V)): NaNs wurden erzeugt

-
## Warning in cov2cor(ans$cov.unscaled): diag(.) had 0 or NA entries; non-
-## finite result is doubtful

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p12a)

+
plot(p12a)

 

-
print(p12a)

+
print(p12a)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 695.4440 220.0685 695.4440 
@@ -609,9 +608,9 @@
 ## $DFOP
 ##          Estimate   Pr(>t)   Lower   Upper
 ## parent_0  100.521 2.74e-10 92.2366 108.805
-## k1          0.124 5.43e-06  0.0959   0.161
-## k2          0.124 6.45e-02  0.0315   0.490
-## g           0.880      NaN      NA      NA
+## k1          0.124 5.74e-06  0.0958   0.161
+## k2          0.124 6.61e-02  0.0319   0.484
+## g           0.877 5.00e-01  0.0000   1.000
 ## sigma       7.048 2.50e-04  4.0349  10.061
 ## 
 ## 
@@ -627,7 +626,7 @@
 

 

 Example on page 12, lower panel

-
p12b <- nafta(NAFTA_SOP_Attachment[["p12b"]])

+
p12b <- nafta(NAFTA_SOP_Attachment[["p12b"]])

 
## Warning in sqrt(diag(covar)): NaNs wurden erzeugt

 
## Warning in qt(alpha/2, rdf): NaNs wurden erzeugt

 
## Warning in qt(1 - alpha/2, rdf): NaNs wurden erzeugt

@@ -638,9 +637,9 @@
 ## finite result is doubtful

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p12b)

+
plot(p12b)

 

-
print(p12b)

+
print(p12b)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 58.90242 19.06353 58.90242 
@@ -667,7 +666,7 @@
 ## parent_0  97.6840    NaN   NaN   NaN
 ## k1         0.0589    NaN    NA    NA
 ## k2         0.0589    NaN    NA    NA
-## g          0.8275    NaN    NA    NA
+## g          0.6902    NaN    NA    NA
 ## sigma      3.4323    NaN   NaN   NaN
 ## 
 ## 
@@ -683,12 +682,16 @@
 

 

 Example on page 13

-
p13 <- nafta(NAFTA_SOP_Attachment[["p13"]])

+
p13 <- nafta(NAFTA_SOP_Attachment[["p13"]])

+
## Warning in sqrt(diag(covar)): NaNs wurden erzeugt

+
## Warning in sqrt(1/diag(V)): NaNs wurden erzeugt

+
## Warning in cov2cor(ans$cov.unscaled): diag(.) had 0 or NA entries; non-
+## finite result is doubtful

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p13)

+
plot(p13)

 

-
print(p13)

+
print(p13)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 174.5971 142.3951 174.5971 
@@ -713,9 +716,9 @@
 ## $DFOP
 ##          Estimate   Pr(>t)    Lower    Upper
 ## parent_0 92.73500 9.25e-15 8.95e+01 9.59e+01
-## k1        0.00258 4.28e-01 1.25e-08 5.31e+02
+## k1        0.00258 4.28e-01 1.70e-08 3.92e+02
 ## k2        0.00258 3.69e-08 2.20e-03 3.03e-03
-## g         0.00442 5.00e-01 0.00e+00 1.00e+00
+## g         0.00442 5.00e-01       NA       NA
 ## sigma     3.41172 1.35e-04 2.02e+00 4.80e+00
 ## 
 ## 
@@ -732,16 +735,16 @@
 

 

 DT50 not observed in the study and DFOP problems in PestDF

-
p14 <- nafta(NAFTA_SOP_Attachment[["p14"]])

+
p14 <- nafta(NAFTA_SOP_Attachment[["p14"]])

 
## Warning in sqrt(diag(covar)): NaNs wurden erzeugt

 
## Warning in sqrt(1/diag(V)): NaNs wurden erzeugt

 
## Warning in cov2cor(ans$cov.unscaled): diag(.) had 0 or NA entries; non-
 ## finite result is doubtful

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p14)

+
plot(p14)

 

-
print(p14)

+
print(p14)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 48.43249 28.67746 27.26248 
@@ -767,7 +770,7 @@
 ##          Estimate   Pr(>t)    Lower    Upper
 ## parent_0 1.00e+02 2.96e-28 99.40280 101.2768
 ## k1       9.53e-03 1.20e-01  0.00638   0.0143
-## k2       5.42e-12 5.00e-01  0.00000      Inf
+## k2       7.29e-12 5.00e-01  0.00000      Inf
 ## g        3.98e-01 2.19e-01  0.30481   0.4998
 ## sigma    1.17e+00 7.68e-06  0.77406   1.5610
 ## 
@@ -776,7 +779,7 @@
 ##          DT50     DT90 DT50_rep
 ## SFO  2.48e+02 8.25e+02 2.48e+02
 ## IORE 4.34e+02 2.22e+04 6.70e+03
-## DFOP 3.41e+10 3.31e+11 1.28e+11
+## DFOP 2.54e+10 2.46e+11 9.51e+10
 ## 
 ## Representative half-life:
 ## [1] 6697.44

@@ -785,7 +788,7 @@
 

 

 N is less than 1 and DFOP fraction parameter is below zero

-
p15a <- nafta(NAFTA_SOP_Attachment[["p15a"]])

+
p15a <- nafta(NAFTA_SOP_Attachment[["p15a"]])

 
## Warning in sqrt(diag(covar)): NaNs wurden erzeugt

 
## Warning in sqrt(diag(covar_notrans)): NaNs wurden erzeugt

 
## Warning in sqrt(1/diag(V)): NaNs wurden erzeugt

@@ -793,9 +796,9 @@
 ## finite result is doubtful

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p15a)

+
plot(p15a)

 

-
print(p15a)

+
print(p15a)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 245.5248 135.0132 245.5248 
@@ -819,10 +822,10 @@
 ## 
 ## $DFOP
 ##          Estimate   Pr(>t)    Lower    Upper
-## parent_0 97.96751 2.85e-13 94.21913 101.7159
-## k1        0.00952 5.68e-02  0.00262   0.0347
-## k2        0.00952 1.52e-04  0.00639   0.0142
-## g         0.22357      NaN       NA       NA
+## parent_0 97.96752 2.85e-13 94.21914 101.7159
+## k1        0.00952 6.80e-02  0.00277   0.0327
+## k2        0.00952 3.82e-06  0.00902   0.0100
+## g         0.17247      NaN       NA       NA
 ## sigma     4.18778 2.50e-04  2.39747   5.9781
 ## 
 ## 
@@ -834,12 +837,16 @@
 ## 
 ## Representative half-life:
 ## [1] 41.33

-
p15b <- nafta(NAFTA_SOP_Attachment[["p15b"]])

+
p15b <- nafta(NAFTA_SOP_Attachment[["p15b"]])

+
## Warning in sqrt(diag(covar)): NaNs wurden erzeugt

+
## Warning in sqrt(1/diag(V)): NaNs wurden erzeugt

+
## Warning in cov2cor(ans$cov.unscaled): diag(.) had 0 or NA entries; non-
+## finite result is doubtful

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The half-life obtained from the IORE model may be used

-
plot(p15b)

+
plot(p15b)

 

-
print(p15b)

+
print(p15b)

 
## Sums of squares:
 ##       SFO      IORE      DFOP 
 ## 106.91629  68.55574 106.91629 
@@ -864,9 +871,9 @@
 ## $DFOP
 ##          Estimate Pr(>t)    Lower    Upper
 ## parent_0 1.01e+02     NA 9.82e+01 1.04e+02
-## k1       4.86e-03     NA 6.49e-04 3.64e-02
-## k2       4.86e-03     NA 3.36e-03 7.03e-03
-## g        1.50e-01     NA 0.00e+00 1.00e+00
+## k1       4.86e-03     NA 6.75e-04 3.49e-02
+## k2       4.86e-03     NA 3.37e-03 6.99e-03
+## g        1.50e-01     NA       NA       NA
 ## sigma    2.76e+00     NA 1.58e+00 3.94e+00
 ## 
 ## 
@@ -883,14 +890,14 @@
 

 

 The DFOP fraction parameter is greater than 1

-
p16 <- nafta(NAFTA_SOP_Attachment[["p16"]])

+
p16 <- nafta(NAFTA_SOP_Attachment[["p16"]])

 
## The SFO model is rejected as S_SFO is equal or higher than the critical value S_c

 
## The representative half-life of the IORE model is longer than the one corresponding

 
## to the terminal degradation rate found with the DFOP model.

 
## The representative half-life obtained from the DFOP model may be used

-
plot(p16)

+
plot(p16)

 

-
print(p16)

+
print(p16)

 
## Sums of squares:
 ##      SFO     IORE     DFOP 
 ## 3831.804 2062.008 1550.980 
@@ -915,7 +922,7 @@
 ## $DFOP
 ##          Estimate   Pr(>t)   Lower  Upper
 ## parent_0  88.5333 7.40e-18 79.9836 97.083
-## k1        18.6315 5.00e-01  0.0000    Inf
+## k1        18.5561 5.00e-01  0.0000    Inf
 ## k2         0.0776 1.41e-05  0.0518  0.116
 ## g          0.4733 1.41e-09  0.3674  0.582
 ## sigma      7.1902 2.11e-08  5.2785  9.102
diff --git a/docs/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png b/docs/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png
index 8ffc5375..396828c3 100644
Binary files a/docs/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png and b/docs/articles/web_only/NAFTA_examples_files/figure-html/p14-1.png differ
diff --git a/docs/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png b/docs/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png
index 48d65737..57bf4014 100644
Binary files a/docs/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png and b/docs/articles/web_only/NAFTA_examples_files/figure-html/p9b-1.png differ
diff --git a/docs/articles/web_only/benchmarks.html b/docs/articles/web_only/benchmarks.html
index 5f6312ce..6df96c14 100644
--- a/docs/articles/web_only/benchmarks.html
+++ b/docs/articles/web_only/benchmarks.html
@@ -88,7 +88,7 @@
       
Benchmark timings for mkin on various systems

                         
Johannes Ranke

             
-            
2019-05-02

+            
2019-05-03

       
       
       
@@ -115,32 +115,32 @@
 }

 
# Parent only
 t1 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), list(FOCUS_2006_C, FOCUS_2006_D)))[["elapsed"]]
-t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), list(FOCUS_2006_C, FOCUS_2006_D), error_model = "tc"))[["elapsed"]]

-
## Warning in mkinfit(models[[model_index]], datasets[[dataset_index]], ...): Optimisation did not converge:
-## false convergence (8)

-
+t2 <- system.time(mmkin_bench(c("SFO", "FOMC", "DFOP", "HS"), list(FOCUS_2006_C, FOCUS_2006_D), error_model = "tc"))[["elapsed"]]
+
+# One metabolite
+SFO_SFO <- mkinmod(
+  parent = mkinsub("SFO", "m1"),
+  m1 = mkinsub("SFO"))

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

-
FOMC_SFO <- mkinmod(
+
+
## Successfully compiled differential equation model from auto-generated C code.

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

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

-
t3 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_2006_D)))[["elapsed"]]

-
## Warning in mkinfit(models[[model_index]], datasets[[dataset_index]], ...):
-## Observations with value of zero were removed from the data

+
t3 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_2006_D)))[["elapsed"]]

 
## Warning in mkinfit(models[[model_index]], datasets[[dataset_index]], ...):
 ## Observations with value of zero were removed from the data
 
+## Warning in mkinfit(models[[model_index]], datasets[[dataset_index]], ...):
+## Observations with value of zero were removed from the data
+
 ## Warning in mkinfit(models[[model_index]], datasets[[dataset_index]], ...):
 ## Observations with value of zero were removed from the data

-
t4 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(subset(FOCUS_2006_D, value != 0)), error_model = "tc"))[["elapsed"]]
-t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_2006_D), error_model = "obs"))[["elapsed"]]

+
t4 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(subset(FOCUS_2006_D, value != 0)), error_model = "tc"))[["elapsed"]]
+t5 <- system.time(mmkin_bench(list(SFO_SFO, FOMC_SFO, DFOP_SFO), list(FOCUS_2006_D), error_model = "obs"))[["elapsed"]]

 
## Warning in mkinfit(models[[model_index]], datasets[[dataset_index]], ...):
 ## Observations with value of zero were removed from the data
 
@@ -149,32 +149,32 @@
 
 ## Warning in mkinfit(models[[model_index]], datasets[[dataset_index]], ...):
 ## Observations with value of zero were removed from the data

-
# Two metabolites, synthetic data
-m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"),
-                           M1 = mkinsub("SFO", "M2"),
-                           M2 = mkinsub("SFO"),
-                           use_of_ff = "max", quiet = TRUE)
-
-m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")),
-                           M1 = mkinsub("SFO"),
-                           M2 = mkinsub("SFO"),
-                           use_of_ff = "max", quiet = TRUE)
-
-SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data
-
-DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data
-
-t6 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a)))["elapsed"]
-t7 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c)))["elapsed"]
-
-t8 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a), error_model = "tc"))["elapsed"]
-t9 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c), error_model = "tc"))["elapsed"]
-
-t10 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a), error_model = "obs"))["elapsed"]
-t11 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c), error_model = "obs"))["elapsed"]
-
-mkin_benchmarks[system_string, paste0("t", 1:11)] <- c(t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11)
-mkin_benchmarks

+
# Two metabolites, synthetic data
+m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"),
+                           M1 = mkinsub("SFO", "M2"),
+                           M2 = mkinsub("SFO"),
+                           use_of_ff = "max", quiet = TRUE)
+
+m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")),
+                           M1 = mkinsub("SFO"),
+                           M2 = mkinsub("SFO"),
+                           use_of_ff = "max", quiet = TRUE)
+
+SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data
+
+DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data
+
+t6 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a)))["elapsed"]
+t7 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c)))["elapsed"]
+
+t8 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a), error_model = "tc"))["elapsed"]
+t9 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c), error_model = "tc"))["elapsed"]
+
+t10 <- system.time(mmkin_bench(list(m_synth_SFO_lin), list(SFO_lin_a), error_model = "obs"))["elapsed"]
+t11 <- system.time(mmkin_bench(list(m_synth_DFOP_par), list(DFOP_par_c), error_model = "obs"))["elapsed"]
+
+mkin_benchmarks[system_string, paste0("t", 1:11)] <- c(t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11)
+mkin_benchmarks

 
##                                                                                                       CPU
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 AMD Ryzen 7 1700 Eight-Core Processor
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 AMD Ryzen 7 1700 Eight-Core Processor
@@ -198,68 +198,68 @@
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 8.184
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 7.064
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 7.296
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 5.850
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 5.864
 ##                                                                         t2
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 11.019
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 22.889
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 12.558
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 21.239
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 20.019
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 23.254
 ##                                                                        t3
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 3.764
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 4.649
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 4.786
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 4.510
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 4.293
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 4.544
 ##                                                                         t4
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 14.347
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 13.789
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2  8.461
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 13.805
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 14.570
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 15.757
 ##                                                                        t5
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 9.495
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 6.395
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 5.675
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 7.386
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 5.772
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 7.870
 ##                                                                        t6
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 2.623
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 2.542
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 2.723
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 2.643
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 2.647
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 2.554
 ##                                                                        t7
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 4.587
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 4.128
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 4.478
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 4.374
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 4.257
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4  4.22
 ##                                                                        t8
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 7.525
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 4.632
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 4.862
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3  7.02
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 4.703
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 6.479
 ##                                                                         t9
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 16.621
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1  8.171
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2  7.618
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 11.124
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4  7.742
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 11.236
 ##                                                                       t10
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 8.576
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1 3.676
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2 3.579
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3 5.388
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 3.421
+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4 4.803
 ##                                                                        t11
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.48.1 31.267
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.1  5.636
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.2  5.574
 ## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.3  7.365
-## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4  5.614

-
save(mkin_benchmarks, file = "~/git/mkin/vignettes/mkin_benchmarks.rda")

+## Linux, AMD Ryzen 7 1700 Eight-Core Processor, mkin version 0.9.49.4  7.688

+
save(mkin_benchmarks, file = "~/git/mkin/vignettes/mkin_benchmarks.rda")

 

   

 
diff --git a/docs/articles/web_only/compiled_models.html b/docs/articles/web_only/compiled_models.html
index 362bc72c..feba6e23 100644
--- a/docs/articles/web_only/compiled_models.html
+++ b/docs/articles/web_only/compiled_models.html
@@ -88,7 +88,7 @@
       
Performance benefit by using compiled model definitions in mkin

                         
Johannes Ranke

             
-            
2019-05-02

+            
2019-05-03

       
       
       
@@ -163,9 +163,9 @@
 ## Warning in mkinfit(SFO_SFO, FOCUS_2006_D, solution_type = "deSolve", quiet
 ## = TRUE): Observations with value of zero were removed from the data

 
##                    test replications elapsed relative user.self sys.self
-## 3     deSolve, compiled            3   3.174    1.000     3.172        0
-## 1 deSolve, not compiled            3  28.777    9.066    28.764        0
-## 2      Eigenvalue based            3   4.386    1.382     4.383        0
+## 3     deSolve, compiled            3   3.075    1.000     3.072    0.000
+## 1 deSolve, not compiled            3  28.192    9.168    28.168    0.008
+## 2      Eigenvalue based            3   4.351    1.415     4.349    0.000
 ##   user.child sys.child
 ## 3          0         0
 ## 1          0         0
@@ -214,8 +214,8 @@
 ## Warning in mkinfit(FOMC_SFO, FOCUS_2006_D, quiet = TRUE): Observations with
 ## value of zero were removed from the data

 
##                    test replications elapsed relative user.self sys.self
-## 2     deSolve, compiled            3   4.561    1.000     4.559        0
-## 1 deSolve, not compiled            3  49.847   10.929    49.823        0
+## 2     deSolve, compiled            3   4.933    1.000     4.930        0
+## 1 deSolve, not compiled            3  52.879   10.719    52.853        0
 ##   user.child sys.child
 ## 2          0         0
 ## 1          0         0

diff --git a/docs/reference/Extract.mmkin.html b/docs/reference/Extract.mmkin.html
index 3b7209b1..652dfee0 100644
--- a/docs/reference/Extract.mmkin.html
+++ b/docs/reference/Extract.mmkin.html
@@ -189,7 +189,7 @@
     fits[["FOMC", "B"]]
   )
#> $par
 #>  parent_0 log_alpha  log_beta     sigma 
-#> 99.666193  2.549850  5.050586  1.890202 
+#> 99.666193  2.549849  5.050586  1.890202 
 #> 
 #> $objective
 #> [1] 28.58291
@@ -202,7 +202,7 @@
 #> 
 #> $evaluations
 #> function gradient 
-#>       25       76 
+#>       25       72 
 #> 
 #> $message
 #> [1] "both X-convergence and relative convergence (5)"
diff --git a/docs/reference/IORE.solution.html b/docs/reference/IORE.solution.html
index c9b79512..c67c2a37 100644
--- a/docs/reference/IORE.solution.html
+++ b/docs/reference/IORE.solution.html
@@ -182,9 +182,9 @@
 
     print(data.frame(fit.fomc$par, fit.iore$par, fit.iore.deS$par,
                      row.names = paste("model par", 1:4)))
#>             fit.fomc.par fit.iore.par fit.iore.deS.par
-#> model par 1  85.87488995    85.874890        85.874890
-#> model par 2   0.05192228    -4.826631        -4.826631
-#> model par 3   0.65096650     1.949403         1.949403
+#> model par 1  85.87489063    85.874890        85.874890
+#> model par 2   0.05192238    -4.826631        -4.826631
+#> model par 3   0.65096665     1.949403         1.949403
 #> model par 4   1.85744396     1.857444         1.857444
    print(rbind(fomc = endpoints(fit.fomc)$distimes, iore = endpoints(fit.iore)$distimes,
                 iore.deS = endpoints(fit.iore)$distimes))
#>              DT50    DT90 DT50back
 #> fomc     1.785233 15.1479 4.559973
diff --git a/docs/reference/NAFTA_SOP_Attachment-1.png b/docs/reference/NAFTA_SOP_Attachment-1.png
index 87068d7c..a40cb689 100644
Binary files a/docs/reference/NAFTA_SOP_Attachment-1.png and b/docs/reference/NAFTA_SOP_Attachment-1.png differ
diff --git a/docs/reference/NAFTA_SOP_Attachment.html b/docs/reference/NAFTA_SOP_Attachment.html
index 48a1495f..76d91c40 100644
--- a/docs/reference/NAFTA_SOP_Attachment.html
+++ b/docs/reference/NAFTA_SOP_Attachment.html
@@ -177,7 +177,7 @@
 #>          Estimate   Pr(>t)   Lower    Upper
 #> parent_0 9.99e+01 1.41e-26 98.8116 101.0810
 #> k1       2.67e-02 5.05e-06  0.0243   0.0295
-#> k2       3.41e-12 5.00e-01  0.0000      Inf
+#> k2       2.86e-12 5.00e-01  0.0000      Inf
 #> g        6.47e-01 3.67e-06  0.6248   0.6677
 #> sigma    1.27e+00 8.91e-06  0.8395   1.6929
 #> 
@@ -186,7 +186,7 @@
 #>      DT50     DT90 DT50_rep
 #> SFO  67.7 2.25e+02 6.77e+01
 #> IORE 58.2 1.07e+03 3.22e+02
-#> DFOP 55.5 3.70e+11 2.03e+11
+#> DFOP 55.5 4.42e+11 2.42e+11
 #> 
 #> Representative half-life:
 #> [1] 321.51
  plot(nafta_att_p5a)
diff --git a/docs/reference/add_err-1.png b/docs/reference/add_err-1.png index 2f60fdf7..61f00d9b 100644 Binary files a/docs/reference/add_err-1.png and b/docs/reference/add_err-1.png differ diff --git a/docs/reference/logLik.mkinfit.html b/docs/reference/logLik.mkinfit.html index 0184d573..f5844a8e 100644 --- a/docs/reference/logLik.mkinfit.html +++ b/docs/reference/logLik.mkinfit.html @@ -180,7 +180,7 @@ The total number of estimated parameters returned with the value f_nw <- mkinfit(sfo_sfo, d_t, quiet = TRUE) # no weighting (weights are unity)
#> Warning: Observations with value of zero were removed from the data
f_obs <- mkinfit(sfo_sfo, d_t, error_model = "obs", quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
f_tc <- mkinfit(sfo_sfo, d_t, error_model = "tc", quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
AIC(f_nw, f_obs, f_tc)
#> df AIC #> f_nw 5 204.4486 #> f_obs 6 205.8727 -#> f_tc 6 141.9656
+#> f_tc 6 148.1802
summary(m)$bpar
#> Estimate se_notrans t value Pr(>t) Lower -#> parent_0 1.057896e+02 1.9023449643 55.610120 3.768361e-16 1.016451e+02 +#> parent_0 1.057896e+02 1.9023449649 55.610120 3.768361e-16 1.016451e+02 #> kmax 6.398190e-02 0.0143201029 4.467978 3.841828e-04 3.929235e-02 -#> k0 1.612775e-04 0.0005866813 0.274898 3.940351e-01 5.846687e-08 -#> r 2.263946e-01 0.1718110718 1.317695 1.061044e-01 4.335843e-02 +#> k0 1.612775e-04 0.0005866813 0.274898 3.940351e-01 5.846685e-08 +#> r 2.263946e-01 0.1718110773 1.317695 1.061044e-01 4.335843e-02 #> sigma 5.332935e+00 0.9145907310 5.830952 4.036926e-05 3.340213e+00 #> Upper #> parent_0 109.9341588 #> kmax 0.1041853 -#> k0 0.4448749 +#> k0 0.4448750 #> r 1.1821121 #> sigma 7.3256566
endpoints(m)$distimes
#> DT50 DT90 DT50_k0 DT50_kmax -#> parent 36.86533 62.41511 4297.853 10.83349
+#> parent 36.86533 62.41511 4297.854 10.83349
#> Compilation argument: -#> /usr/lib/R/bin/R CMD SHLIB file26cc3b169ccf.c 2> file26cc3b169ccf.c.err.txt +#> /usr/lib/R/bin/R CMD SHLIB file2b304be466d0.c 2> file2b304be466d0.c.err.txt #> Program source: #> 1: #include <R.h> #> 2: diff --git a/docs/reference/mkinparplot-1.png b/docs/reference/mkinparplot-1.png index 798722a2..65918d53 100644 Binary files a/docs/reference/mkinparplot-1.png and b/docs/reference/mkinparplot-1.png differ diff --git a/docs/reference/mkinparplot.html b/docs/reference/mkinparplot.html index c90acf4a..8cfb9374 100644 --- a/docs/reference/mkinparplot.html +++ b/docs/reference/mkinparplot.html @@ -154,8 +154,7 @@ 
model <- mkinmod(
   T245 = mkinsub("SFO", to = c("phenol"), sink = FALSE),
   phenol = mkinsub("SFO", to = c("anisole")),
-  anisole = mkinsub("SFO"), use_of_ff = "max")
#> Successfully compiled differential equation model from auto-generated C code.
fit <- mkinfit(model, subset(mccall81_245T, soil == "Commerce"), quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
#> Warning: Optimisation did not converge:
-#> false convergence (8)
mkinparplot(fit)
+ anisole = mkinsub("SFO"), use_of_ff = "max")
#> Successfully compiled differential equation model from auto-generated C code.
fit <- mkinfit(model, subset(mccall81_245T, soil == "Commerce"), quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
mkinparplot(fit)
#> time parent m1 #> 201 20 4.978707 27.46227
#> User System verstrichen -#> 0.002 0.000 0.002
system.time( +#> 0.002 0.000 0.001
system.time( print(mkinpredict(SFO_SFO, c(k_parent_m1 = 0.05, k_parent_sink = 0.1, k_m1_sink = 0.01), c(parent = 100, m1 = 0), seq(0, 20, by = 0.1), solution_type = "deSolve", use_compiled = FALSE)[201,]))
#> time parent m1 #> 201 20 4.978707 27.46227
#> User System verstrichen -#> 0.021 0.000 0.021
+#> 0.022 0.000 0.021
# Predict from a fitted model - f <- mkinfit(SFO_SFO, FOCUS_2006_C)
#> Sum of squared residuals at call 1: 552.5739 + f <- mkinfit(SFO_SFO, FOCUS_2006_C)
#> Ordinary least squares optimisation
#> Sum of squared residuals at call 1: 552.5739 #> Sum of squared residuals at call 3: 552.5739 #> Sum of squared residuals at call 4: 552.5739 #> Sum of squared residuals at call 6: 279.9345 @@ -350,6 +350,7 @@ #> Sum of squared residuals at call 12: 200.3629 #> Sum of squared residuals at call 13: 200.3629 #> Sum of squared residuals at call 18: 197.9039 +#> Sum of squared residuals at call 23: 197.9039 #> Sum of squared residuals at call 25: 196.6754 #> Sum of squared residuals at call 27: 196.6754 #> Sum of squared residuals at call 32: 196.5742 @@ -366,15 +367,16 @@ #> Sum of squared residuals at call 58: 196.5334 #> Sum of squared residuals at call 59: 196.5334 #> Sum of squared residuals at call 65: 196.5334 +#> Sum of squared residuals at call 73: 196.5334 #> Negative log-likelihood at call 75: 26.64668 #> Negative log-likelihood at call 103: 26.64668 #> Optimisation successfully terminated.
head(mkinpredict(f))
#> time parent m1 #> 1 0.0 82.49216 0.000000 -#> 2 0.1 80.00563 1.179956 -#> 3 0.2 77.59404 2.312583 -#> 4 0.3 75.25515 3.399424 -#> 5 0.4 72.98675 4.441974 -#> 6 0.5 70.78673 5.441685
+#> 2 0.1 80.00563 1.179955 +#> 3 0.2 77.59404 2.312580 +#> 4 0.3 75.25515 3.399419 +#> 5 0.4 72.98675 4.441969 +#> 6 0.5 70.78673 5.441679
#> User System verstrichen -#> 0.047 0.041 7.417
time_1
#> User System verstrichen -#> 21.251 0.000 21.263
+#> 0.042 0.028 5.095
time_1
#> User System verstrichen +#> 19.29 0.00 19.30
endpoints(fits.0[["SFO_lin", 2]])
#> $ff #> parent_M1 parent_sink M1_M2 M1_sink -#> 0.7340481 0.2659519 0.7505683 0.2494317 +#> 0.7340481 0.2659519 0.7505684 0.2494316 #> #> $SFORB #> logical(0) @@ -206,8 +206,8 @@ #> $distimes #> DT50 DT90 #> parent 0.8777689 2.915885 -#> M1 2.3257442 7.725955 -#> M2 33.7200941 112.015728 +#> M1 2.3257449 7.725957 +#> M2 33.7200958 112.015734 #>
# plot.mkinfit handles rows or columns of mmkin result objects plot(fits.0[1, ])
plot(fits.0[1, ], obs_var = c("M1", "M2"))
plot(fits.0[, 1])
# Use double brackets to extract a single mkinfit object, which will be plotted diff --git a/docs/reference/plot.mkinfit-1.png b/docs/reference/plot.mkinfit-1.png index 115792db..e7b3bac7 100644 Binary files a/docs/reference/plot.mkinfit-1.png and b/docs/reference/plot.mkinfit-1.png differ diff --git a/docs/reference/plot.mmkin-3.png b/docs/reference/plot.mmkin-3.png index 412c7244..7f32afe2 100644 Binary files a/docs/reference/plot.mmkin-3.png and b/docs/reference/plot.mmkin-3.png differ diff --git a/docs/reference/schaefer07_complex_case.html b/docs/reference/schaefer07_complex_case.html index f645e6a7..09f6d1e2 100644 --- a/docs/reference/schaefer07_complex_case.html +++ b/docs/reference/schaefer07_complex_case.html @@ -167,7 +167,7 @@ A2 = list(type = "SFO"), use_of_ff = "max")
#> Successfully compiled differential equation model from auto-generated C code.
fit <- mkinfit(model, data, quiet = TRUE) plot(fit)
endpoints(fit)
#> $ff #> parent_A1 parent_B1 parent_C1 parent_sink A1_A2 A1_sink -#> 0.3809619 0.1954667 0.4235714 0.0000000 0.4479603 0.5520397 +#> 0.3809619 0.1954667 0.4235714 0.0000000 0.4479620 0.5520380 #> #> $SFORB #> logical(0) @@ -175,10 +175,10 @@ #> $distimes #> DT50 DT90 #> parent 13.95078 46.34350 -#> A1 49.75343 165.27733 -#> B1 37.26907 123.80518 -#> C1 11.23131 37.30959 -#> A2 28.50644 94.69635 +#> A1 49.75344 165.27737 +#> B1 37.26907 123.80517 +#> C1 11.23131 37.30960 +#> A2 28.50627 94.69577 #>
# Compare with the results obtained in the original publication print(schaefer07_complex_results)
#> compound parameter KinGUI ModelMaker deviation #> 1 parent degradation rate 0.0496 0.0506 2.0 diff --git a/docs/reference/summary.mkinfit.html b/docs/reference/summary.mkinfit.html index 89a8b9a6..351d7a7e 100644 --- a/docs/reference/summary.mkinfit.html +++ b/docs/reference/summary.mkinfit.html @@ -211,15 +211,15 @@

Examples

  summary(mkinfit(mkinmod(parent = mkinsub("SFO")), FOCUS_2006_A, quiet = TRUE))
#> mkin version used for fitting:    0.9.49.4 
 #> R version used for fitting:       3.6.0 
-#> Date of fit:     Thu May  2 18:48:59 2019 
-#> Date of summary: Thu May  2 18:48:59 2019 
+#> Date of fit:     Fri May  3 19:08:31 2019 
+#> Date of summary: Fri May  3 19:08:31 2019 
 #> 
 #> Equations:
 #> d_parent/dt = - k_parent_sink * parent
 #> 
 #> Model predictions using solution type analytical 
 #> 
-#> Fitted using 131 model solutions performed in 0.274 s
+#> Fitted using 131 model solutions performed in 0.27 s
 #> 
 #> Error model:
 #> Constant variance 
@@ -247,9 +247,9 @@
 #> 
 #> Parameter correlation:
 #>                    parent_0 log_k_parent_sink     sigma
-#> parent_0          1.000e+00         5.428e-01 1.642e-07
-#> log_k_parent_sink 5.428e-01         1.000e+00 2.507e-07
-#> sigma             1.642e-07         2.507e-07 1.000e+00
+#> parent_0          1.000e+00         5.428e-01 1.648e-07
+#> log_k_parent_sink 5.428e-01         1.000e+00 2.513e-07
+#> sigma             1.648e-07         2.513e-07 1.000e+00
 #> 
 #> Backtransformed parameters:
 #> Confidence intervals for internally transformed parameters are asymmetric.
diff --git a/docs/reference/test_data_from_UBA_2014.html b/docs/reference/test_data_from_UBA_2014.html
index 688baf32..bc988340 100644
--- a/docs/reference/test_data_from_UBA_2014.html
+++ b/docs/reference/test_data_from_UBA_2014.html
@@ -160,14 +160,14 @@
   m_ws <- mkinmod(parent_w = mkinsub("SFO", "parent_s"),
                   parent_s = mkinsub("SFO", "parent_w"))
#> Successfully compiled differential equation model from auto-generated C code.
  f_river <- mkinfit(m_ws, test_data_from_UBA_2014[[1]]$data, quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
  plot_sep(f_river)

   summary(f_river)$bpar
#>                         Estimate se_notrans      t value       Pr(>t)
-#> parent_w_0          9.598567e+01 2.12352039 4.520120e+01 9.476357e-18
-#> k_parent_w_sink     3.603743e-01 0.03149366 1.144276e+01 4.128096e-09
-#> k_parent_w_parent_s 6.031371e-02 0.01603582 3.761185e+00 9.436293e-04
-#> k_parent_s_sink     7.560341e-11 0.09483761 7.971881e-10 5.000000e-01
-#> k_parent_s_parent_w 7.419672e-02 0.10738374 6.909493e-01 2.500756e-01
-#> sigma               2.982879e+00 0.50546582 5.901247e+00 1.454824e-05
+#> parent_w_0          9.598567e+01 2.12351789 4.520126e+01 9.476190e-18
+#> k_parent_w_sink     3.603743e-01 0.03149282 1.144306e+01 4.126593e-09
+#> k_parent_w_parent_s 6.031371e-02 0.01603582 3.761186e+00 9.436275e-04
+#> k_parent_s_sink     5.108964e-11 0.09482736 5.387647e-10 5.000000e-01
+#> k_parent_s_parent_w 7.419672e-02 0.10737376 6.910135e-01 2.500560e-01
+#> sigma               2.982879e+00 0.50545649 5.901356e+00 1.454535e-05
 #>                           Lower       Upper
-#> parent_w_0          91.48420501 100.4871438
+#> parent_w_0          91.48420503 100.4871438
 #> k_parent_w_sink      0.30668904   0.4234571
 #> k_parent_w_parent_s  0.03423904   0.1062455
 #> k_parent_s_sink      0.00000000         Inf
@@ -184,25 +184,25 @@
                     M3 = mkinsub("SFO"),
                     use_of_ff = "max")
#> Successfully compiled differential equation model from auto-generated C code.

   f_soil <- mkinfit(m_soil, test_data_from_UBA_2014[[3]]$data, quiet = TRUE)
#> Warning: Observations with value of zero were removed from the data
  plot_sep(f_soil, lpos = c("topright", "topright", "topright", "bottomright"))
  summary(f_soil)$bpar
#>                   Estimate  se_notrans    t value       Pr(>t)        Lower
-#> parent_0       76.55425583 0.859186612 89.1008482 1.113866e-26 74.755959748
-#> k_parent        0.12081956 0.004601921 26.2541551 1.077372e-16  0.111561582
-#> k_M1            0.84258650 0.806231456  1.0450926 1.545475e-01  0.113839804
-#> k_M2            0.04210878 0.017083049  2.4649452 1.170195e-02  0.018013807
+#> parent_0       76.55425584 0.859186619 89.1008474 1.113866e-26 74.755959756
+#> k_parent        0.12081956 0.004601922 26.2541544 1.077373e-16  0.111561582
+#> k_M1            0.84258649 0.806231419  1.0450926 1.545475e-01  0.113839803
+#> k_M2            0.04210878 0.017083049  2.4649453 1.170195e-02  0.018013807
 #> k_M3            0.01122919 0.007245890  1.5497322 6.885127e-02  0.002909463
-#> f_parent_to_M1  0.32240199 0.240803564  1.3388589 9.820820e-02           NA
-#> f_parent_to_M2  0.16099854 0.033691991  4.7785403 6.531225e-05           NA
-#> f_M1_to_M3      0.27921500 0.269443517  1.0362654 1.565440e-01  0.022992933
-#> f_M2_to_M3      0.55641333 0.595125466  0.9349513 1.807725e-01  0.008003317
+#> f_parent_to_M1  0.32240199 0.240803555  1.3388589 9.820820e-02           NA
+#> f_parent_to_M2  0.16099854 0.033691991  4.7785403 6.531224e-05           NA
+#> f_M1_to_M3      0.27921501 0.269443514  1.0362655 1.565440e-01  0.022992937
+#> f_M2_to_M3      0.55641331 0.595125445  0.9349513 1.807725e-01  0.008003317
 #> sigma           1.14005399 0.149696423  7.6157731 1.727024e-07  0.826735778
 #>                      Upper
 #> parent_0       78.35255192
 #> k_parent        0.13084582
-#> k_M1            6.23641283
+#> k_M1            6.23641265
 #> k_M2            0.09843279
 #> k_M3            0.04333950
 #> f_parent_to_M1          NA
 #> f_parent_to_M2          NA
-#> f_M1_to_M3      0.86443084
+#> f_M1_to_M3      0.86443083
 #> f_M2_to_M3      0.99489847
 #> sigma           1.45337221
  mkinerrmin(f_soil)
#>             err.min n.optim df
 #> All data 0.09649963       9 20
diff --git a/test.log b/test.log
index b1e9b9f3..bd6c6d0e 100644
--- a/test.log
+++ b/test.log
@@ -2,26 +2,26 @@ Loading mkin
 Testing mkin
 ✔ |  OK F W S | Context
 
⠏ |   0       | Export dataset for reading into CAKE
✔ |   1       | Export dataset for reading into CAKE
-
⠏ |   0       | Error model fitting
⠋ |   1       | Error model fitting
⠹ |   3       | Error model fitting
⠸ |   4       | Error model fitting
⠼ |   5       | Error model fitting
⠴ |   6       | Error model fitting
⠧ |   8       | Error model fitting
⠏ |  10       | Error model fitting
⠋ |  11       | Error model fitting
✔ |  12       | Error model fitting [147.8 s]
-
⠏ |   0       | Calculation of FOCUS chi2 error levels
⠋ |   1       | Calculation of FOCUS chi2 error levels
⠹ |   3       | Calculation of FOCUS chi2 error levels
✔ |   3       | Calculation of FOCUS chi2 error levels [2.3 s]
-
⠏ |   0       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠋ |   1       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠙ |   2       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠸ |   4       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠇ |   9       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
✔ |  13       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.7 s]
+
⠏ |   0       | Error model fitting
⠋ |   1       | Error model fitting
⠹ |   3       | Error model fitting
⠸ |   4       | Error model fitting
⠼ |   5       | Error model fitting
⠴ |   6       | Error model fitting
⠧ |   8       | Error model fitting
⠏ |  10       | Error model fitting
⠋ |  11       | Error model fitting
✔ |  12       | Error model fitting [219.8 s]
+
⠏ |   0       | Calculation of FOCUS chi2 error levels
⠋ |   1       | Calculation of FOCUS chi2 error levels
⠹ |   3       | Calculation of FOCUS chi2 error levels
✔ |   3       | Calculation of FOCUS chi2 error levels [2.4 s]
+
⠏ |   0       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠋ |   1       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠙ |   2       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠸ |   4       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
⠇ |   9       | Results for FOCUS D established in expertise for UBA (Ranke 2014)
✔ |  13       | Results for FOCUS D established in expertise for UBA (Ranke 2014) [4.0 s]
 
⠏ |   0       | Test fitting the decline of metabolites from their maximum
⠋ |   1       | Test fitting the decline of metabolites from their maximum
⠹ |   3       | Test fitting the decline of metabolites from their maximum
⠼ |   5       | Test fitting the decline of metabolites from their maximum
✔ |   6       | Test fitting the decline of metabolites from their maximum [0.9 s]
 
⠏ |   0       | Fitting the logistic model
⠋ |   1       | Fitting the logistic model
✔ |   1       | Fitting the logistic model [0.9 s]
 
⠏ |   0       | Test dataset class mkinds used in gmkin
✔ |   1       | Test dataset class mkinds used in gmkin
-
⠏ |   0       | Special cases of mkinfit calls
⠋ |   1       | Special cases of mkinfit calls
⠇ |   9       | Special cases of mkinfit calls
⠏ |  10       | Special cases of mkinfit calls
⠋ |  11       | Special cases of mkinfit calls
⠙ |  12       | Special cases of mkinfit calls
✔ |  12       | Special cases of mkinfit calls [2.7 s]
+
⠏ |   0       | Special cases of mkinfit calls
⠋ |   1       | Special cases of mkinfit calls
⠇ |   9       | Special cases of mkinfit calls
⠏ |  10       | Special cases of mkinfit calls
⠋ |  11       | Special cases of mkinfit calls
⠙ |  12       | Special cases of mkinfit calls
✔ |  12       | Special cases of mkinfit calls [2.9 s]
 
⠏ |   0       | mkinmod model generation and printing
⠇ |   9       | mkinmod model generation and printing
✔ |   9       | mkinmod model generation and printing [0.2 s]
 
⠏ |   0       | Model predictions with mkinpredict
⠋ |   1       | Model predictions with mkinpredict
✔ |   3       | Model predictions with mkinpredict [0.3 s]
-
⠏ |   0       | Evaluations according to 2015 NAFTA guidance
⠙ |   2       | Evaluations according to 2015 NAFTA guidance
⠇ |   9       | Evaluations according to 2015 NAFTA guidance
⠏ |  10       | Evaluations according to 2015 NAFTA guidance
⠴ |  16       | Evaluations according to 2015 NAFTA guidance
✔ |  16       | Evaluations according to 2015 NAFTA guidance [3.9 s]
-
⠏ |   0       | Fitting of parent only models
⠋ |   1       | Fitting of parent only models
⠙ |   2       | Fitting of parent only models
⠹ |   3       | Fitting of parent only models
⠸ |   4       | Fitting of parent only models
⠼ |   5       | Fitting of parent only models
⠴ |   6       | Fitting of parent only models
⠦ |   7       | Fitting of parent only models
⠧ |   8       | Fitting of parent only models
⠇ |   9       | Fitting of parent only models
⠏ |  10       | Fitting of parent only models
⠋ |  11       | Fitting of parent only models
⠙ |  12       | Fitting of parent only models
⠹ |  13       | Fitting of parent only models
⠴ |  16       | Fitting of parent only models
⠧ |  18       | Fitting of parent only models
⠏ |  20       | Fitting of parent only models
✔ |  21       | Fitting of parent only models [40.1 s]
+
⠏ |   0       | Evaluations according to 2015 NAFTA guidance
⠙ |   2       | Evaluations according to 2015 NAFTA guidance
⠇ |   9       | Evaluations according to 2015 NAFTA guidance
⠏ |  10       | Evaluations according to 2015 NAFTA guidance
⠴ |  16       | Evaluations according to 2015 NAFTA guidance
✔ |  16       | Evaluations according to 2015 NAFTA guidance [4.0 s]
+
⠏ |   0       | Fitting of parent only models
⠋ |   1       | Fitting of parent only models
⠙ |   2       | Fitting of parent only models
⠹ |   3       | Fitting of parent only models
⠸ |   4       | Fitting of parent only models
⠼ |   5       | Fitting of parent only models
⠴ |   6       | Fitting of parent only models
⠦ |   7       | Fitting of parent only models
⠧ |   8       | Fitting of parent only models
⠇ |   9       | Fitting of parent only models
⠏ |  10       | Fitting of parent only models
⠋ |  11       | Fitting of parent only models
⠙ |  12       | Fitting of parent only models
⠹ |  13       | Fitting of parent only models
⠴ |  16       | Fitting of parent only models
⠧ |  18       | Fitting of parent only models
⠏ |  20       | Fitting of parent only models
✔ |  21       | Fitting of parent only models [40.8 s]
 
⠏ |   0       | Calculation of maximum time weighted average concentrations (TWAs)
⠋ |   1       | Calculation of maximum time weighted average concentrations (TWAs)
⠙ |   2       | Calculation of maximum time weighted average concentrations (TWAs)
⠹ |   3       | Calculation of maximum time weighted average concentrations (TWAs)
⠸ |   4       | Calculation of maximum time weighted average concentrations (TWAs)
✔ |   4       | Calculation of maximum time weighted average concentrations (TWAs) [2.2 s]
 
⠏ |   0       | Summary
✔ |   1       | Summary
 
⠏ |   0       | Plotting
⠹ |   3       | Plotting
✔ |   4       | Plotting [0.3 s]
 
⠏ |   0       | AIC calculation
✔ |   2       | AIC calculation
-
⠏ |   0       | Complex test case from Schaefer et al. (2007) Piacenza paper
⠋ |   1       | Complex test case from Schaefer et al. (2007) Piacenza paper
✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [5.3 s]
+
⠏ |   0       | Complex test case from Schaefer et al. (2007) Piacenza paper
⠋ |   1       | Complex test case from Schaefer et al. (2007) Piacenza paper
✔ |   2       | Complex test case from Schaefer et al. (2007) Piacenza paper [5.4 s]
 
⠏ |   0       | Results for synthetic data established in expertise for UBA (Ranke 2014)
⠋ |   1       | Results for synthetic data established in expertise for UBA (Ranke 2014)
⠹ |   3       | Results for synthetic data established in expertise for UBA (Ranke 2014)
✔ |   4       | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.1 s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 223.4 s
+Duration: 297.0 s
 
 OK:       115
 Failed:   0
diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf
index 07488d5e..02c7893f 100644
--- a/tests/testthat/FOCUS_2006_D.csf
+++ b/tests/testthat/FOCUS_2006_D.csf
@@ -5,7 +5,7 @@ Description:
 MeasurementUnits: % AR
 TimeUnits: days
 Comments: Created using mkin::CAKE_export
-Date: 2019-05-02
+Date: 2019-05-03
 Optimiser: IRLS
 
 [Data]
diff --git a/vignettes/mkin_benchmarks.rda b/vignettes/mkin_benchmarks.rda
index 61f03620..985cdac6 100644
Binary files a/vignettes/mkin_benchmarks.rda and b/vignettes/mkin_benchmarks.rda differ
-- 
cgit v1.2.1