From ce73c044b949154e3bc3e715b9b79e1360b3f794 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 1 Nov 2019 15:34:28 +0100 Subject: Make the 'quadratic' the default for 'confint' Also the documentation was improved here and there --- DESCRIPTION | 2 +- NEWS.md | 4 + R/AIC.mmkin.R | 8 +- R/confint.mkinfit.R | 64 ++++--- R/parms.mkinfit.R | 6 +- R/residuals.mkinfit.R | 2 +- R/update.mkinfit.R | 8 +- docs/404.html | 2 +- docs/articles/index.html | 2 +- docs/authors.html | 2 +- docs/index.html | 2 +- docs/news/index.html | 11 +- docs/reference/AIC.mmkin.html | 17 +- docs/reference/confint.mkinfit.html | 109 ++++++----- docs/reference/index.html | 2 +- docs/reference/parms.html | 19 +- docs/reference/residuals.mkinfit.html | 4 +- docs/reference/update.mkinfit-1.png | Bin 0 -> 28678 bytes docs/reference/update.mkinfit-2.png | Bin 0 -> 28541 bytes docs/reference/update.mkinfit.html | 350 +--------------------------------- man/AIC.mmkin.Rd | 8 +- man/confint.mkinfit.Rd | 66 ++++--- man/parms.Rd | 5 + man/residuals.mkinfit.Rd | 2 +- man/update.mkinfit.Rd | 8 +- 25 files changed, 232 insertions(+), 471 deletions(-) create mode 100644 docs/reference/update.mkinfit-1.png create mode 100644 docs/reference/update.mkinfit-2.png diff --git a/DESCRIPTION b/DESCRIPTION index d204bdd3..b0bcc39d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data -Version: 0.9.49.7 +Version: 0.9.49.8 Date: 2019-11-01 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de", diff --git a/NEWS.md b/NEWS.md index 78423592..b9d58a27 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# mkin 0.9.49.8 (unreleased) + +- 'confint.mkinfit': Make the quadratic approximation the default, as the likelihood profiling takes too much time for any with more than three parameters + # mkin 0.9.49.7 (2019-11-01) - Fix a bug introduced in 0.9.49.6 that occurred if the direct optimisation yielded a higher likelihood than the three-step optimisation in the d_3 algorithm, which caused the fitted parameters of the three-step optimisation to be returned instead of the parameters of the direct optimisation diff --git a/R/AIC.mmkin.R b/R/AIC.mmkin.R index 7d405c4a..f1a66998 100644 --- a/R/AIC.mmkin.R +++ b/R/AIC.mmkin.R @@ -17,15 +17,21 @@ #' f <- mmkin(c("SFO", "FOMC", "DFOP"), #' list("FOCUS A" = FOCUS_2006_A, #' "FOCUS C" = FOCUS_2006_C), cores = 1, quiet = TRUE) -#' AIC(f[1, "FOCUS A"]) # We get a single number for a single fit +#' # We get a warning because the FOMC model does not converge for the +#' # FOCUS A dataset, as it is well described by SFO +#' +#' AIC(f["SFO", "FOCUS A"]) # We get a single number for a single fit +#' AIC(f[["SFO", "FOCUS A"]]) # or when extracting an mkinfit object #' #' # For FOCUS A, the models fit almost equally well, so the higher the number #' # of parameters, the higher (worse) the AIC #' AIC(f[, "FOCUS A"]) #' AIC(f[, "FOCUS A"], k = 0) # If we do not penalize additional parameters, we get nearly the same +#' BIC(f[, "FOCUS A"]) # Comparing the BIC gives a very similar picture #' #' # For FOCUS C, the more complex models fit better #' AIC(f[, "FOCUS C"]) +#' BIC(f[, "FOCUS C"]) #' } #' #' @export diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R index 5e1703d6..07614306 100644 --- a/R/confint.mkinfit.R +++ b/R/confint.mkinfit.R @@ -1,8 +1,13 @@ #' Confidence intervals for parameters of mkinfit objects #' -#' The default method 'profile' is based on the profile likelihood for each -#' parameter. The method uses two nested optimisations. The speed of the method -#' could likely be improved by using the method of Venzon and Moolgavkar (1988). +#' The default method 'quadratic' is based on the quadratic approximation of +#' the curvature of the likelihood function at the maximum likelihood parameter +#' estimates. +#' The alternative method 'profile' is based on the profile likelihood for each +#' parameter. The method uses two nested optimisations and can take a very long +#' time, even if parallelized by specifying 'cores' on unixoid platforms. The +#' speed of the method could likely be improved by using the method of Venzon +#' and Moolgavkar (1988). #' #' @param object An \code{\link{mkinfit}} object #' @param parm A vector of names of the parameters which are to be given @@ -12,11 +17,11 @@ #' @param cutoff Possibility to specify an alternative cutoff for the difference #' in the log-likelihoods at the confidence boundary. Specifying an explicit #' cutoff value overrides arguments 'level' and 'alpha' -#' @param method The 'profile' method searches the parameter space for the +#' @param method The 'quadratic' method approximates the likelihood function at +#' the optimised parameters using the second term of the Taylor expansion, +#' using a second derivative (hessian) contained in the object. +#' The 'profile' method searches the parameter space for the #' cutoff of the confidence intervals by means of a likelihood ratio test. -#' The 'quadratic' method approximates the likelihood function at the -#' optimised parameters using the second term of the Taylor expansion, using -#' a second derivative (hessian) contained in the object. #' @param transformed If the quadratic approximation is used, should it be #' applied to the likelihood based on the transformed parameters? #' @param backtransform If we approximate the likelihood in terms of the @@ -46,19 +51,26 @@ #' \dontrun{ #' confint(f, method = "profile") #' +#' # Set the number of cores for the profiling method for further examples +#' if (identical(Sys.getenv("NOT_CRAN"), "true")) { +#' n_cores <- parallel::detectCores() - 1 +#' } else { +#' n_cores <- 1 +#' } +#' if (Sys.getenv("TRAVIS") != "") n_cores = 1 +#' if (Sys.info()["sysname"] == "Windows") n_cores = 1 +#' #' SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), quiet = TRUE) #' SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), #' use_of_ff = "max", quiet = TRUE) #' f_d_1 <- mkinfit(SFO_SFO, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -#' system.time(ci_profile <- confint(f_d_1, cores = 1, quiet = TRUE)) -#' # The following does not save much time, as parent_0 takes up most of the time -#' # system.time(ci_profile <- confint(f_d_1, cores = 5)) -#' # system.time(ci_profile <- confint(f_d_1, -#' # c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 1)) -#' # If we exclude parent_0 (the confidence of which is often of minor interest), we get a nice -#' # performance improvement from about 30 seconds to about 12 seconds -#' # system.time(ci_profile_no_parent_0 <- confint(f_d_1, -#' # c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 4)) +#' system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE)) +#' # Using more cores does not save much time here, as parent_0 takes up most of the time +#' # If we additionally exclude parent_0 (the confidence of which is often of +#' # minor interest), we get a nice performance improvement from about 50 +#' # seconds to about 12 seconds if we use at least four cores +#' system.time(ci_profile_no_parent_0 <- confint(f_d_1, method = "profile", +#' c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores)) #' ci_profile #' ci_quadratic_transformed <- confint(f_d_1, method = "quadratic") #' ci_quadratic_transformed @@ -70,21 +82,14 @@ #' # interval based on the untransformed fit for k_m1_sink #' rel_diffs_transformed <- abs((ci_quadratic_transformed - ci_profile)/ci_profile) #' rel_diffs_untransformed <- abs((ci_quadratic_untransformed - ci_profile)/ci_profile) -#' rel_diffs_transformed -#' rel_diffs_untransformed +#' rel_diffs_transformed < rel_diffs_untransformed +#' signif(rel_diffs_transformed, 3) +#' signif(rel_diffs_untransformed, 3) #' -#' # Set the number of cores for further examples -#' if (identical(Sys.getenv("NOT_CRAN"), "true")) { -#' n_cores <- parallel::detectCores() - 1 -#' } else { -#' n_cores <- 1 -#' } -#' if (Sys.getenv("TRAVIS") != "") n_cores = 1 -#' if (Sys.info()["sysname"] == "Windows") n_cores = 1 #' #' # Investigate a case with formation fractions #' f_d_2 <- mkinfit(SFO_SFO.ff, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -#' ci_profile_ff <- confint(f_d_2, cores = n_cores) +#' ci_profile_ff <- confint(f_d_2, method = "profile", cores = n_cores) #' ci_profile_ff #' ci_quadratic_transformed_ff <- confint(f_d_2, method = "quadratic") #' ci_quadratic_transformed_ff @@ -94,8 +99,9 @@ #' rel_diffs_untransformed_ff <- abs((ci_quadratic_untransformed_ff - ci_profile_ff)/ci_profile_ff) #' # While the confidence interval for the parent rate constant is closer to #' # the profile based interval when using the internal parameter -#' # transformation, the intervals for the other parameters are 'better +#' # transformation, the interval for the metabolite rate constant is 'better #' # without internal parameter transformation. +#' rel_diffs_transformed_ff < rel_diffs_untransformed_ff #' rel_diffs_transformed_ff #' rel_diffs_untransformed_ff #' @@ -114,7 +120,7 @@ #' @export confint.mkinfit <- function(object, parm, level = 0.95, alpha = 1 - level, cutoff, - method = c("profile", "quadratic"), + method = c("quadratic", "profile"), transformed = TRUE, backtransform = TRUE, cores = round(detectCores()/2), quiet = FALSE, ...) { diff --git a/R/parms.mkinfit.R b/R/parms.mkinfit.R index 0628cb92..281d06de 100644 --- a/R/parms.mkinfit.R +++ b/R/parms.mkinfit.R @@ -3,7 +3,7 @@ #' This function always returns degradation model parameters as well as error #' model parameters, in order to avoid working with a fitted model without #' considering the error structure that was assumed for the fit. -#' +#' #' @param object A fitted model object #' @param \dots Not used #' @return A numeric vector of fitted model parameters @@ -16,6 +16,10 @@ parms <- function(object, ...) #' @param transformed Should the parameters be returned #' as used internally during the optimisation? #' @rdname parms +#' @examples +#' fit <- mkinfit("SFO", FOCUS_2006_C) +#' parms(fit) +#' parms(fit, transformed = TRUE) #' @export parms.mkinfit <- function(object, transformed = FALSE, ...) { diff --git a/R/residuals.mkinfit.R b/R/residuals.mkinfit.R index 96bcf01c..2a508fc3 100644 --- a/R/residuals.mkinfit.R +++ b/R/residuals.mkinfit.R @@ -1,6 +1,6 @@ #' Extract residuals from an mkinfit model #' -#' @param object An \code{\link{mkinfit}} object +#' @param object A \code{\link{mkinfit}} object #' @param standardized Should the residuals be standardized by dividing by the #' standard deviation obtained from the fitted error model? #' @param \dots Not used diff --git a/R/update.mkinfit.R b/R/update.mkinfit.R index 2f0814e0..dde7f810 100644 --- a/R/update.mkinfit.R +++ b/R/update.mkinfit.R @@ -12,8 +12,12 @@ #' @param evaluate Should the call be evaluated or returned as a call #' @examples #' \dontrun{ -#' fit <- mkinfit("DFOP", subset(FOCUS_2006_D, value != 0), quiet = TRUE) -#' update(fit, error_model = "tc") +#' fit <- mkinfit("SFO", subset(FOCUS_2006_D, value != 0), quiet = TRUE) +#' parms(fit) +#' plot_err(fit) +#' fit_2 <- update(fit, error_model = "tc") +#' parms(fit_2) +#' plot_err(fit_2) #' } #' @export update.mkinfit <- function(object, ..., evaluate = TRUE) diff --git a/docs/404.html b/docs/404.html index f10cadb9..2a71e496 100644 --- a/docs/404.html +++ b/docs/404.html @@ -67,7 +67,7 @@ mkin - 0.9.49.7 + 0.9.49.8 diff --git a/docs/articles/index.html b/docs/articles/index.html index 76d7bc69..ac5f5df1 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -67,7 +67,7 @@ mkin - 0.9.49.7 + 0.9.49.8 diff --git a/docs/authors.html b/docs/authors.html index cb970788..fb50f268 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -67,7 +67,7 @@ mkin - 0.9.49.7 + 0.9.49.8 diff --git a/docs/index.html b/docs/index.html index 5a618e57..231ef39f 100644 --- a/docs/index.html +++ b/docs/index.html @@ -38,7 +38,7 @@ mkin - 0.9.49.7 + 0.9.49.8 diff --git a/docs/news/index.html b/docs/news/index.html index 5d1fc00f..baddfc96 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -67,7 +67,7 @@ mkin - 0.9.49.7 + 0.9.49.8 @@ -129,6 +129,14 @@ +
+

+mkin 0.9.49.8 (unreleased) Unreleased +

+ +

mkin 0.9.49.7 (2019-11-01) Unreleased @@ -735,6 +743,7 @@

Contents

    +
  • 0.9.49.8
  • 0.9.49.7
  • 0.9.49.6
  • 0.9.49.5
    • diff --git a/docs/reference/AIC.mmkin.html b/docs/reference/AIC.mmkin.html index 7ca12302..6b782456 100644 --- a/docs/reference/AIC.mmkin.html +++ b/docs/reference/AIC.mmkin.html @@ -70,7 +70,7 @@ same dataset." /> mkin - 0.9.49.7 + 0.9.49.8
@@ -173,7 +173,10 @@ column.

f <- mmkin(c("SFO", "FOMC", "DFOP"), list("FOCUS A" = FOCUS_2006_A, "FOCUS C" = FOCUS_2006_C), cores = 1, quiet = TRUE)

#> Warning: Optimisation did not converge: -#> false convergence (8)
AIC(f[1, "FOCUS A"]) # We get a single number for a single fit
#> [1] 55.28197
+#> false convergence (8)
# We get a warning because the FOMC model does not converge for the + # FOCUS A dataset, as it is well described by SFO + + AIC(f["SFO", "FOCUS A"]) # We get a single number for a single fit
#> [1] 55.28197
AIC(f[["SFO", "FOCUS A"]]) # or when extracting an mkinfit object
#> [1] 55.28197
# For FOCUS A, the models fit almost equally well, so the higher the number # of parameters, the higher (worse) the AIC AIC(f[, "FOCUS A"])
#> df AIC @@ -182,12 +185,18 @@ column.

#> DFOP 5 59.28197
AIC(f[, "FOCUS A"], k = 0) # If we do not penalize additional parameters, we get nearly the same
#> df AIC #> SFO 3 49.28197 #> FOMC 4 49.28202 -#> DFOP 5 49.28197
+#> DFOP 5 49.28197
BIC(f[, "FOCUS A"]) # Comparing the BIC gives a very similar picture
#> df BIC +#> SFO 3 55.52030 +#> FOMC 4 57.59979 +#> DFOP 5 59.67918
# For FOCUS C, the more complex models fit better AIC(f[, "FOCUS C"])
#> df AIC #> SFO 3 59.29336 #> FOMC 4 44.68652 -#> DFOP 5 29.02372
+#> DFOP 5 29.02372
BIC(f[, "FOCUS C"])
#> df BIC +#> SFO 3 59.88504 +#> FOMC 4 45.47542 +#> DFOP 5 30.00984
diff --git a/docs/reference/confint.mkinfit.html b/docs/reference/confint.mkinfit.html index 54696ff5..6015fbf6 100644 --- a/docs/reference/confint.mkinfit.html +++ b/docs/reference/confint.mkinfit.html @@ -36,9 +36,14 @@ - + @@ -71,7 +76,7 @@ could likely be improved by using the method of Venzon and Moolgavkar (1988)." / mkin - 0.9.49.6 + 0.9.49.8 @@ -135,14 +140,19 @@ could likely be improved by using the method of Venzon and Moolgavkar (1988)." /
-

The default method 'profile' is based on the profile likelihood for each -parameter. The method uses two nested optimisations. The speed of the method -could likely be improved by using the method of Venzon and Moolgavkar (1988).

+

The default method 'quadratic' is based on the quadratic approximation of +the curvature of the likelihood function at the maximum likelihood parameter +estimates. +The alternative method 'profile' is based on the profile likelihood for each +parameter. The method uses two nested optimisations and can take a very long +time, even if parallelized by specifying 'cores' on unixoid platforms. The +speed of the method could likely be improved by using the method of Venzon +and Moolgavkar (1988).

# S3 method for mkinfit
 confint(object, parm, level = 0.95, alpha = 1 -
-  level, cutoff, method = c("profile", "quadratic"),
+  level, cutoff, method = c("quadratic", "profile"),
   transformed = TRUE, backtransform = TRUE,
   cores = round(detectCores()/2), quiet = FALSE, ...)
@@ -174,11 +184,11 @@ cutoff value overrides arguments 'level' and 'alpha'

method -

The 'profile' method searches the parameter space for the -cutoff of the confidence intervals by means of a likelihood ratio test. -The 'quadratic' method approximates the likelihood function at the -optimised parameters using the second term of the Taylor expansion, using -a second derivative (hessian) contained in the object.

+

The 'quadratic' method approximates the likelihood function at +the optimised parameters using the second term of the Taylor expansion, +using a second derivative (hessian) contained in the object. +The 'profile' method searches the parameter space for the +cutoff of the confidence intervals by means of a likelihood ratio test.

transformed @@ -231,20 +241,27 @@ machines, cores > 1 is not supported.

#> parent_0 73.0641834 92.1392181 #> k_parent_sink 0.2170293 0.4235348 #> sigma 3.1307772 8.0628314
+# Set the number of cores for the profiling method for further examples +if (identical(Sys.getenv("NOT_CRAN"), "true")) { + n_cores <- parallel::detectCores() - 1 +} else { + n_cores <- 1 +} +if (Sys.getenv("TRAVIS") != "") n_cores = 1 +if (Sys.info()["sysname"] == "Windows") n_cores = 1 + SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), quiet = TRUE) SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) f_d_1 <- mkinfit(SFO_SFO, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -system.time(ci_profile <- confint(f_d_1, cores = 1, quiet = TRUE))
#> User System verstrichen -#> 51.646 0.000 51.673
# The following does not save much time, as parent_0 takes up most of the time -# system.time(ci_profile <- confint(f_d_1, cores = 5)) -# system.time(ci_profile <- confint(f_d_1, -# c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 1)) -# If we exclude parent_0 (the confidence of which is often of minor interest), we get a nice -# performance improvement from about 30 seconds to about 12 seconds -# system.time(ci_profile_no_parent_0 <- confint(f_d_1, -# c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 4)) -ci_profile
#> 2.5% 97.5% +system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE))
#> User System verstrichen +#> 51.297 0.000 51.328
# Using more cores does not save much time here, as parent_0 takes up most of the time +# If we additionally exclude parent_0 (the confidence of which is often of +# minor interest), we get a nice performance improvement from about 50 +# seconds to about 12 seconds if we use at least four cores +system.time(ci_profile_no_parent_0 <- confint(f_d_1, method = "profile", + c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores))
#> Profiling the likelihood
#> User System verstrichen +#> 0.006 0.003 11.435
ci_profile
#> 2.5% 97.5% #> parent_0 96.456003650 1.027703e+02 #> k_parent_sink 0.040762501 5.549764e-02 #> k_parent_m1 0.046786482 5.500879e-02 @@ -267,29 +284,26 @@ machines, cores > 1 is not supported.

# interval based on the untransformed fit for k_m1_sink rel_diffs_transformed <- abs((ci_quadratic_transformed - ci_profile)/ci_profile) rel_diffs_untransformed <- abs((ci_quadratic_untransformed - ci_profile)/ci_profile) -rel_diffs_transformed
#> 2.5% 97.5% -#> parent_0 0.0005407854 0.0002218012 -#> k_parent_sink 0.0066452394 0.0083795930 -#> k_parent_m1 0.0001833903 0.0020092090 -#> k_m1_sink 0.0307278240 0.0290580487 -#> sigma 0.0550252516 0.0327066836
rel_diffs_untransformed
#> 2.5% 97.5% -#> parent_0 0.0005407854 0.0002218011 -#> k_parent_sink 0.0067996407 0.0025717594 -#> k_parent_m1 0.0037382781 0.0011843074 -#> k_m1_sink 0.0146745610 0.0025299672 -#> sigma 0.0550252516 0.0327066836
-# Set the number of cores for further examples -if (identical(Sys.getenv("NOT_CRAN"), "true")) { - n_cores <- parallel::detectCores() - 1 -} else { - n_cores <- 1 -} -if (Sys.getenv("TRAVIS") != "") n_cores = 1 -if (Sys.info()["sysname"] == "Windows") n_cores = 1 +rel_diffs_transformed < rel_diffs_untransformed
#> 2.5% 97.5% +#> parent_0 FALSE FALSE +#> k_parent_sink TRUE FALSE +#> k_parent_m1 TRUE FALSE +#> k_m1_sink FALSE FALSE +#> sigma FALSE FALSE
signif(rel_diffs_transformed, 3)
#> 2.5% 97.5% +#> parent_0 0.000541 0.000222 +#> k_parent_sink 0.006650 0.008380 +#> k_parent_m1 0.000183 0.002010 +#> k_m1_sink 0.030700 0.029100 +#> sigma 0.055000 0.032700
signif(rel_diffs_untransformed, 3)
#> 2.5% 97.5% +#> parent_0 0.000541 0.000222 +#> k_parent_sink 0.006800 0.002570 +#> k_parent_m1 0.003740 0.001180 +#> k_m1_sink 0.014700 0.002530 +#> sigma 0.055000 0.032700
# Investigate a case with formation fractions f_d_2 <- mkinfit(SFO_SFO.ff, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -ci_profile_ff <- confint(f_d_2, cores = n_cores)
#> Profiling the likelihood
ci_profile_ff
#> 2.5% 97.5% +ci_profile_ff <- confint(f_d_2, method = "profile", cores = n_cores)
#> Profiling the likelihood
ci_profile_ff
#> 2.5% 97.5% #> parent_0 96.456003650 1.027703e+02 #> k_parent 0.090911032 1.071578e-01 #> k_m1 0.003892605 6.702778e-03 @@ -310,9 +324,14 @@ machines, cores > 1 is not supported.

rel_diffs_untransformed_ff <- abs((ci_quadratic_untransformed_ff - ci_profile_ff)/ci_profile_ff) # While the confidence interval for the parent rate constant is closer to # the profile based interval when using the internal parameter -# transformation, the intervals for the other parameters are 'better +# transformation, the interval for the metabolite rate constant is 'better # without internal parameter transformation. -rel_diffs_transformed_ff
#> 2.5% 97.5% +rel_diffs_transformed_ff < rel_diffs_untransformed_ff
#> 2.5% 97.5% +#> parent_0 TRUE TRUE +#> k_parent TRUE TRUE +#> k_m1 FALSE FALSE +#> f_parent_to_m1 TRUE FALSE +#> sigma FALSE TRUE
rel_diffs_transformed_ff
#> 2.5% 97.5% #> parent_0 0.0005408012 0.0002217857 #> k_parent 0.0009596303 0.0009003981 #> k_m1 0.0307277425 0.0290579163 diff --git a/docs/reference/index.html b/docs/reference/index.html index 4279f6de..1c9975e0 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -67,7 +67,7 @@ mkin - 0.9.49.7 + 0.9.49.8
diff --git a/docs/reference/parms.html b/docs/reference/parms.html index 8ab26240..e3c52c53 100644 --- a/docs/reference/parms.html +++ b/docs/reference/parms.html @@ -71,7 +71,7 @@ considering the error structure that was assumed for the fit." /> mkin - 0.9.49.6 + 0.9.49.8 @@ -167,12 +167,29 @@ as used internally during the optimisation?

A numeric vector of fitted model parameters

+

Examples

+ 
fit <- mkinfit("SFO", FOCUS_2006_C)
#> Ordinary least squares optimisation
#> Sum of squared residuals at call 1: 2388.077
+#> Sum of squared residuals at call 3: 2388.077
+#> Sum of squared residuals at call 4: 247.1962
+#> Sum of squared residuals at call 7: 200.6791
+#> Sum of squared residuals at call 10: 197.7231
+#> Sum of squared residuals at call 11: 197.0872
+#> Sum of squared residuals at call 14: 196.535
+#> Sum of squared residuals at call 15: 196.535
+#> Sum of squared residuals at call 16: 196.535
+#> Sum of squared residuals at call 17: 196.5334
+#> Sum of squared residuals at call 20: 196.5334
+#> Sum of squared residuals at call 25: 196.5334
+#> Negative log-likelihood at call 31: 26.64668
#> Optimisation successfully terminated.
parms(fit)
#>      parent_0 k_parent_sink         sigma 
+#>    82.4921598     0.3060633     4.6730124 
parms(fit, transformed = TRUE)
#>          parent_0 log_k_parent_sink             sigma 
+#>         82.492160         -1.183963          4.673012
diff --git a/docs/reference/residuals.mkinfit.html b/docs/reference/residuals.mkinfit.html index 7dd25012..96b02908 100644 --- a/docs/reference/residuals.mkinfit.html +++ b/docs/reference/residuals.mkinfit.html @@ -69,7 +69,7 @@ mkin - 0.9.49.6 + 0.9.49.7 @@ -144,7 +144,7 @@ object -

An mkinfit object

+

A mkinfit object

standardized diff --git a/docs/reference/update.mkinfit-1.png b/docs/reference/update.mkinfit-1.png new file mode 100644 index 00000000..4cbefa56 Binary files /dev/null and b/docs/reference/update.mkinfit-1.png differ diff --git a/docs/reference/update.mkinfit-2.png b/docs/reference/update.mkinfit-2.png new file mode 100644 index 00000000..f432f6f8 Binary files /dev/null and b/docs/reference/update.mkinfit-2.png differ diff --git a/docs/reference/update.mkinfit.html b/docs/reference/update.mkinfit.html index b9f6f9a9..bb93a61e 100644 --- a/docs/reference/update.mkinfit.html +++ b/docs/reference/update.mkinfit.html @@ -72,7 +72,7 @@ override these starting values." /> mkin - 0.9.49.6 + 0.9.49.8 @@ -167,349 +167,11 @@ remove arguments given in the original call

Examples

# \dontrun{
-fit <- mkinfit("DFOP", subset(FOCUS_2006_D, value != 0), quiet = TRUE)
-update(fit, error_model = "tc")
#> $par
-#>     parent_0       log_k1       log_k2        g_ilr    sigma_low     rsd_high 
-#> 100.85489822   8.84468217  -2.29693632 -14.95263998   0.00375222   0.06763435 
-#> 
-#> $objective
-#> [1] 19.40656
-#> 
-#> $convergence
-#> [1] 0
-#> 
-#> $iterations
-#> [1] 120
-#> 
-#> $evaluations
-#> function gradient 
-#>      144      847 
-#> 
-#> $message
-#> [1] "relative convergence (4)"
-#> 
-#> $logLik
-#> [1] -19.40656
-#> 
-#> $d_3_message
-#>                                                            threestep 
-#> "Three-step fitting yielded a higher likelihood than direct fitting" 
-#> 
-#> $hessian
-#>                parent_0        log_k1        log_k2         g_ilr     sigma_low
-#> parent_0   3.662473e-01 -2.914408e-16 -7.241561e+01 -3.021629e-08  1.923504e+01
-#> log_k1    -2.914408e-16  0.000000e+00 -2.215935e-13 -7.291307e-25  3.416474e-15
-#> log_k2    -7.241561e+01 -2.215935e-13  3.127457e+04  6.766544e-06 -1.495826e+04
-#> g_ilr     -3.021629e-08 -7.291307e-25  6.766544e-06  3.122099e-09 -1.797429e-06
-#> sigma_low  1.923504e+01  3.416474e-15 -1.495826e+04 -1.797429e-06  7.759299e+04
-#> rsd_high   3.902119e+00 -1.801019e-16 -1.685343e+02 -3.750713e-07  3.984179e+03
-#>                rsd_high
-#> parent_0   3.902119e+00
-#> log_k1    -1.801019e-16
-#> log_k2    -1.685343e+02
-#> g_ilr     -3.750713e-07
-#> sigma_low  3.984179e+03
-#> rsd_high   7.188991e+03
-#> 
-#> $hessian_notrans
-#>                parent_0            k1            k2             g     sigma_low
-#> parent_0   3.662473e-01 -3.714445e-19 -7.201669e+02 -3.261485e+01  1.923504e+01
-#> k1        -3.714445e-19  0.000000e+00 -4.258512e-15  7.218123e-21  4.355854e-18
-#> k2        -7.201669e+02 -4.258512e-15  3.092510e+06  7.263235e+04 -1.456870e+05
-#> g         -3.261485e+01  7.218123e-21  7.263235e+04  3.291750e+03 -1.939948e+03
-#> sigma_low  1.923504e+01  4.355854e-18 -1.456870e+05 -1.939948e+03  7.759299e+04
-#> rsd_high   3.902119e+00 -2.259812e-19 -1.779680e+03 -4.048658e+02  3.984179e+03
-#>                rsd_high
-#> parent_0   3.902119e+00
-#> k1        -2.259812e-19
-#> k2        -1.779680e+03
-#> g         -4.048658e+02
-#> sigma_low  3.984179e+03
-#> rsd_high   7.188991e+03
-#> 
-#> $call
-#> mkinfit(mkinmod = "DFOP", observed = subset(FOCUS_2006_D, value != 
-#>     0), parms.ini = c(k1 = 0.699298911979803, k2 = 0.0899931270871125, 
-#> g = 0.0923391681138686), state.ini = c(parent = 101.948852047129), 
-#>     quiet = TRUE, error_model = "tc")
-#> 
-#> $error_model_algorithm
-#> [1] "d_3"
-#> 
-#> $solution_type
-#> [1] "analytical"
-#> 
-#> $transform_rates
-#> [1] TRUE
-#> 
-#> $transform_fractions
-#> [1] TRUE
-#> 
-#> $reweight.tol
-#> [1] 1e-08
-#> 
-#> $reweight.max.iter
-#> [1] 10
-#> 
-#> $control
-#> $control$eval.max
-#> [1] 300
-#> 
-#> $control$iter.max
-#> [1] 200
-#> 
-#> 
-#> $calls
-#> [1] 3105
-#> 
-#> $time
-#>        User      System verstrichen 
-#>      10.251       0.000      10.257 
-#> 
-#> $mkinmod
-#> <mkinmod> model generated with
-#> Use of formation fractions $use_of_ff: min 
-#> Specification $spec:
-#> $parent
-#> $type: DFOP; $sink: TRUE
-#> Differential equations:
-#> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
-#>            time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time)))
-#>            * parent
-#> 
-#> $observed
-#>      name time  value
-#> 1  parent    0  99.46
-#> 2  parent    0 102.04
-#> 3  parent    1  93.50
-#> 4  parent    1  92.50
-#> 5  parent    3  63.23
-#> 6  parent    3  68.99
-#> 7  parent    7  52.32
-#> 8  parent    7  55.13
-#> 9  parent   14  27.27
-#> 10 parent   14  26.64
-#> 11 parent   21  11.50
-#> 12 parent   21  11.64
-#> 13 parent   35   2.85
-#> 14 parent   35   2.91
-#> 15 parent   50   0.69
-#> 16 parent   50   0.63
-#> 17 parent   75   0.05
-#> 18 parent   75   0.06
-#> 
-#> $obs_vars
-#> [1] "parent"
-#> 
-#> $predicted
-#>       name       time        value
-#> 1   parent  0.0000000 100.85489822
-#> 2   parent  0.7575758  93.45650191
-#> 3   parent  1.0000000  91.20560523
-#> 4   parent  1.5151515  86.60082860
-#> 5   parent  2.2727273  80.24806580
-#> 6   parent  3.0000000  74.58828105
-#> 7   parent  3.0303030  74.36132158
-#> 8   parent  3.7878788  68.90641029
-#> 9   parent  4.5454545  63.85165403
-#> 10  parent  5.3030303  59.16769870
-#> 11  parent  6.0606061  54.82734352
-#> 12  parent  6.8181818  50.80538306
-#> 13  parent  7.0000000  49.88485755
-#> 14  parent  7.5757576  47.07846089
-#> 15  parent  8.3333333  43.62493394
-#> 16  parent  9.0909091  40.42474681
-#> 17  parent  9.8484848  37.45931528
-#> 18  parent 10.6060606  34.71141842
-#> 19  parent 11.3636364  32.16509858
-#> 20  parent 12.1212121  29.80556871
-#> 21  parent 12.8787879  27.61912649
-#> 22  parent 13.6363636  25.59307475
-#> 23  parent 14.0000000  24.67405211
-#> 24  parent 14.3939394  23.71564776
-#> 25  parent 15.1515152  21.97594287
-#> 26  parent 15.9090909  20.36385722
-#> 27  parent 16.6666667  18.87002909
-#> 28  parent 17.4242424  17.48578345
-#> 29  parent 18.1818182  16.20308170
-#> 30  parent 18.9393939  15.01447489
-#> 31  parent 19.6969697  13.91306051
-#> 32  parent 20.4545455  12.89244241
-#> 33  parent 21.0000000  12.20428157
-#> 34  parent 21.2121212  11.94669362
-#> 35  parent 21.9696970  11.07032198
-#> 36  parent 22.7272727  10.25823818
-#> 37  parent 23.4848485   9.50572628
-#> 38  parent 24.2424242   8.80841627
-#> 39  parent 25.0000000   8.16225872
-#> 40  parent 25.7575758   7.56350125
-#> 41  parent 26.5151515   7.00866672
-#> 42  parent 27.2727273   6.49453311
-#> 43  parent 28.0303030   6.01811471
-#> 44  parent 28.7878788   5.57664485
-#> 45  parent 29.5454545   5.16755983
-#> 46  parent 30.3030303   4.78848399
-#> 47  parent 31.0606061   4.43721595
-#> 48  parent 31.8181818   4.11171583
-#> 49  parent 32.5757576   3.81009336
-#> 50  parent 33.3333333   3.53059697
-#> 51  parent 34.0909091   3.27160354
-#> 52  parent 34.8484848   3.03160906
-#> 53  parent 35.0000000   2.98576554
-#> 54  parent 35.6060606   2.80921981
-#> 55  parent 36.3636364   2.60314433
-#> 56  parent 37.1212121   2.41218590
-#> 57  parent 37.8787879   2.23523557
-#> 58  parent 38.6363636   2.07126576
-#> 59  parent 39.3939394   1.91932426
-#> 60  parent 40.1515152   1.77852870
-#> 61  parent 40.9090909   1.64806147
-#> 62  parent 41.6666667   1.52716489
-#> 63  parent 42.4242424   1.41513691
-#> 64  parent 43.1818182   1.31132694
-#> 65  parent 43.9393939   1.21513214
-#> 66  parent 44.6969697   1.12599389
-#> 67  parent 45.4545455   1.04339454
-#> 68  parent 46.2121212   0.96685442
-#> 69  parent 46.9696970   0.89592904
-#> 70  parent 47.7272727   0.83020652
-#> 71  parent 48.4848485   0.76930520
-#> 72  parent 49.2424242   0.71287140
-#> 73  parent 50.0000000   0.66057741
-#> 74  parent 50.7575758   0.61211954
-#> 75  parent 51.5151515   0.56721639
-#> 76  parent 52.2727273   0.52560719
-#> 77  parent 53.0303030   0.48705031
-#> 78  parent 53.7878788   0.45132184
-#> 79  parent 54.5454545   0.41821430
-#> 80  parent 55.3030303   0.38753542
-#> 81  parent 56.0606061   0.35910705
-#> 82  parent 56.8181818   0.33276409
-#> 83  parent 57.5757576   0.30835357
-#> 84  parent 58.3333333   0.28573373
-#> 85  parent 59.0909091   0.26477320
-#> 86  parent 59.8484848   0.24535028
-#> 87  parent 60.6060606   0.22735216
-#> 88  parent 61.3636364   0.21067432
-#> 89  parent 62.1212121   0.19521992
-#> 90  parent 62.8787879   0.18089921
-#> 91  parent 63.6363636   0.16762901
-#> 92  parent 64.3939394   0.15533228
-#> 93  parent 65.1515152   0.14393759
-#> 94  parent 65.9090909   0.13337879
-#> 95  parent 66.6666667   0.12359454
-#> 96  parent 67.4242424   0.11452804
-#> 97  parent 68.1818182   0.10612662
-#> 98  parent 68.9393939   0.09834151
-#> 99  parent 69.6969697   0.09112749
-#> 100 parent 70.4545455   0.08444266
-#> 101 parent 71.2121212   0.07824822
-#> 102 parent 71.9696970   0.07250818
-#> 103 parent 72.7272727   0.06718920
-#> 104 parent 73.4848485   0.06226042
-#> 105 parent 74.2424242   0.05769319
-#> 106 parent 75.0000000   0.05346100
-#> 
-#> $rss
-#> function (P) 
-#> cost_function(P, OLS = TRUE, update_data = FALSE)
-#> <bytecode: 0x555558e48258>
-#> <environment: 0x55555caa08e0>
-#> 
-#> $ll
-#> function (P, fixed_degparms = FALSE, fixed_errparms = FALSE) 
-#> {
-#>     -cost_function(P, trans = FALSE, fixed_degparms = fixed_degparms, 
-#>         fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE)
-#> }
-#> <bytecode: 0x555558e47e30>
-#> <environment: 0x55555caa08e0>
-#> 
-#> $start
-#>                  value   type
-#> parent_0  101.94885205  state
-#> k1          0.69929891 deparm
-#> k2          0.08999313 deparm
-#> g           0.09233917 deparm
-#> sigma_low   0.10000000  error
-#> rsd_high    0.10000000  error
-#> 
-#> $start_transformed
-#>                value lower upper
-#> parent_0  101.948852  -Inf   Inf
-#> log_k1     -0.357677  -Inf   Inf
-#> log_k2     -2.408022  -Inf   Inf
-#> g_ilr      -1.616024  -Inf   Inf
-#> sigma_low   0.100000     0   Inf
-#> rsd_high    0.100000     0   Inf
-#> 
-#> $fixed
-#> [1] value type 
-#> <0 Zeilen> (oder row.names mit Länge 0)
-#> 
-#> $data
-#>    time variable observed   predicted     residual
-#> 1     0   parent    99.46 100.8548982  -1.39489822
-#> 2     0   parent   102.04 100.8548982   1.18510178
-#> 3     1   parent    93.50  91.2056052   2.29439477
-#> 4     1   parent    92.50  91.2056052   1.29439477
-#> 5     3   parent    63.23  74.5882810 -11.35828105
-#> 6     3   parent    68.99  74.5882810  -5.59828105
-#> 7     7   parent    52.32  49.8848576   2.43514245
-#> 8     7   parent    55.13  49.8848576   5.24514245
-#> 9    14   parent    27.27  24.6740521   2.59594789
-#> 10   14   parent    26.64  24.6740521   1.96594789
-#> 11   21   parent    11.50  12.2042816  -0.70428157
-#> 12   21   parent    11.64  12.2042816  -0.56428157
-#> 13   35   parent     2.85   2.9857655  -0.13576554
-#> 14   35   parent     2.91   2.9857655  -0.07576554
-#> 15   50   parent     0.69   0.6605774   0.02942259
-#> 16   50   parent     0.63   0.6605774  -0.03057741
-#> 17   75   parent     0.05   0.0534610  -0.00346100
-#> 18   75   parent     0.06   0.0534610   0.00653900
-#> 
-#> $atol
-#> [1] 1e-08
-#> 
-#> $rtol
-#> [1] 1e-10
-#> 
-#> $err_mod
-#> [1] "tc"
-#> 
-#> $bparms.optim
-#>     parent_0           k1           k2            g 
-#> 1.008549e+02 6.937399e+03 1.005665e-01 6.551046e-10 
-#> 
-#> $bparms.fixed
-#> numeric(0)
-#> 
-#> $bparms.ode
-#>           k1           k2            g 
-#> 6.937399e+03 1.005665e-01 6.551046e-10 
-#> 
-#> $bparms.state
-#>   parent 
-#> 100.8549 
-#> 
-#> $errparms
-#>  sigma_low   rsd_high 
-#> 0.00375222 0.06763435 
-#> 
-#> $df.residual
-#> [1] 12
-#> 
-#> $date
-#> [1] "Fri Nov  1 10:10:42 2019"
-#> 
-#> $version
-#> [1] "0.9.49.6"
-#> 
-#> $Rversion
-#> [1] "3.6.1"
-#> 
-#> attr(,"class")
-#> [1] "mkinfit" "modFit" 
# }
+fit <- mkinfit("SFO", subset(FOCUS_2006_D, value != 0), quiet = TRUE)
+parms(fit)
#>      parent_0 k_parent_sink         sigma 
+#>   99.44423886    0.09793574    3.39632469 
plot_err(fit)
fit_2 <- update(fit, error_model = "tc")
+parms(fit_2)
#>      parent_0 k_parent_sink     sigma_low      rsd_high 
+#>  1.008549e+02  1.005665e-01  3.752222e-03  6.763434e-02 
plot_err(fit_2)
# }