diff options
| author | Johannes Ranke <jranke@uni-bremen.de> | 2019-11-01 15:34:28 +0100 |
|---|---|---|
| committer | Johannes Ranke <jranke@uni-bremen.de> | 2019-11-01 15:34:28 +0100 |
| commit | ce73c044b949154e3bc3e715b9b79e1360b3f794 (patch) | |
| tree | 28f477a3142f1efa463a8d8569f924b9064e8637 | |
| parent | e7c65ee913d4a84da0957d2ebb89abfbc444de56 (diff) | |
Make the 'quadratic' the default for 'confint'
Also the documentation was improved here and there
| -rw-r--r-- | DESCRIPTION | 2 | ||||
| -rw-r--r-- | NEWS.md | 4 | ||||
| -rw-r--r-- | R/AIC.mmkin.R | 8 | ||||
| -rw-r--r-- | R/confint.mkinfit.R | 64 | ||||
| -rw-r--r-- | R/parms.mkinfit.R | 6 | ||||
| -rw-r--r-- | R/residuals.mkinfit.R | 2 | ||||
| -rw-r--r-- | R/update.mkinfit.R | 8 | ||||
| -rw-r--r-- | docs/404.html | 2 | ||||
| -rw-r--r-- | docs/articles/index.html | 2 | ||||
| -rw-r--r-- | docs/authors.html | 2 | ||||
| -rw-r--r-- | docs/index.html | 2 | ||||
| -rw-r--r-- | docs/news/index.html | 11 | ||||
| -rw-r--r-- | docs/reference/AIC.mmkin.html | 17 | ||||
| -rw-r--r-- | docs/reference/confint.mkinfit.html | 109 | ||||
| -rw-r--r-- | docs/reference/index.html | 2 | ||||
| -rw-r--r-- | docs/reference/parms.html | 19 | ||||
| -rw-r--r-- | docs/reference/residuals.mkinfit.html | 4 | ||||
| -rw-r--r-- | docs/reference/update.mkinfit-1.png | bin | 0 -> 28678 bytes | |||
| -rw-r--r-- | docs/reference/update.mkinfit-2.png | bin | 0 -> 28541 bytes | |||
| -rw-r--r-- | docs/reference/update.mkinfit.html | 350 | ||||
| -rw-r--r-- | man/AIC.mmkin.Rd | 8 | ||||
| -rw-r--r-- | man/confint.mkinfit.Rd | 66 | ||||
| -rw-r--r-- | man/parms.Rd | 5 | ||||
| -rw-r--r-- | man/residuals.mkinfit.Rd | 2 | ||||
| -rw-r--r-- | man/update.mkinfit.Rd | 8 |
25 files changed, 232 insertions, 471 deletions
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", @@ -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 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="https://pkgdown.jrwb.de/mkin/index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> 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 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> 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 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> 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 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> 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 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> @@ -129,6 +129,14 @@ </div> + <div id="mkin-0-9-49-8-unreleased" class="section level1"> +<h1 class="page-header"> +<a href="#mkin-0-9-49-8-unreleased" class="anchor"></a>mkin 0.9.49.8 (unreleased)<small> Unreleased </small> +</h1> +<ul> +<li>‘confint.mkinfit’: Make the quadratic approximation the default, as the likelihood profiling takes too much time for any with more than three parameters</li> +</ul> +</div> <div id="mkin-0-9-49-7-2019-11-01" class="section level1"> <h1 class="page-header"> <a href="#mkin-0-9-49-7-2019-11-01" class="anchor"></a>mkin 0.9.49.7 (2019-11-01)<small> Unreleased </small> @@ -735,6 +743,7 @@ <div id="tocnav"> <h2>Contents</h2> <ul class="nav nav-pills nav-stacked"> + <li><a href="#mkin-0-9-49-8-unreleased">0.9.49.8</a></li> <li><a href="#mkin-0-9-49-7-2019-11-01">0.9.49.7</a></li> <li><a href="#mkin-0-9-49-6-2019-10-31">0.9.49.6</a></li> <li><a href="#mkin-0-9-49-5-2019-07-04">0.9.49.5</a></li> 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." /> </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> @@ -173,7 +173,10 @@ column.</p></td> <span class='no'>f</span> <span class='kw'><-</span> <span class='fu'><a href='mmkin.html'>mmkin</a></span>(<span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span>(<span class='st'>"SFO"</span>, <span class='st'>"FOMC"</span>, <span class='st'>"DFOP"</span>), <span class='fu'><a href='https://rdrr.io/r/base/list.html'>list</a></span>(<span class='st'>"FOCUS A"</span> <span class='kw'>=</span> <span class='no'>FOCUS_2006_A</span>, <span class='st'>"FOCUS C"</span> <span class='kw'>=</span> <span class='no'>FOCUS_2006_C</span>), <span class='kw'>cores</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)</div><div class='output co'>#> <span class='warning'>Warning: Optimisation did not converge:</span> -#> <span class='warning'>false convergence (8)</span></div><div class='input'> <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>AIC</a></span>(<span class='no'>f</span>[<span class='fl'>1</span>, <span class='st'>"FOCUS A"</span>]) <span class='co'># We get a single number for a single fit</span></div><div class='output co'>#> [1] 55.28197</div><div class='input'> +#> <span class='warning'>false convergence (8)</span></div><div class='input'> <span class='co'># We get a warning because the FOMC model does not converge for the</span> + <span class='co'># FOCUS A dataset, as it is well described by SFO</span> + + <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>AIC</a></span>(<span class='no'>f</span>[<span class='st'>"SFO"</span>, <span class='st'>"FOCUS A"</span>]) <span class='co'># We get a single number for a single fit</span></div><div class='output co'>#> [1] 55.28197</div><div class='input'> <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>AIC</a></span>(<span class='no'>f</span><span class='kw'>[[</span><span class='st'>"SFO"</span>, <span class='st'>"FOCUS A"</span>]]) <span class='co'># or when extracting an mkinfit object</span></div><div class='output co'>#> [1] 55.28197</div><div class='input'> <span class='co'># For FOCUS A, the models fit almost equally well, so the higher the number</span> <span class='co'># of parameters, the higher (worse) the AIC</span> <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>AIC</a></span>(<span class='no'>f</span>[, <span class='st'>"FOCUS A"</span>])</div><div class='output co'>#> df AIC @@ -182,12 +185,18 @@ column.</p></td> #> DFOP 5 59.28197</div><div class='input'> <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>AIC</a></span>(<span class='no'>f</span>[, <span class='st'>"FOCUS A"</span>], <span class='kw'>k</span> <span class='kw'>=</span> <span class='fl'>0</span>) <span class='co'># If we do not penalize additional parameters, we get nearly the same</span></div><div class='output co'>#> df AIC #> SFO 3 49.28197 #> FOMC 4 49.28202 -#> DFOP 5 49.28197</div><div class='input'> +#> DFOP 5 49.28197</div><div class='input'> <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>BIC</a></span>(<span class='no'>f</span>[, <span class='st'>"FOCUS A"</span>]) <span class='co'># Comparing the BIC gives a very similar picture</span></div><div class='output co'>#> df BIC +#> SFO 3 55.52030 +#> FOMC 4 57.59979 +#> DFOP 5 59.67918</div><div class='input'> <span class='co'># For FOCUS C, the more complex models fit better</span> <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>AIC</a></span>(<span class='no'>f</span>[, <span class='st'>"FOCUS C"</span>])</div><div class='output co'>#> df AIC #> SFO 3 59.29336 #> FOMC 4 44.68652 -#> DFOP 5 29.02372</div><div class='input'> +#> DFOP 5 29.02372</div><div class='input'> <span class='fu'><a href='https://rdrr.io/r/stats/AIC.html'>BIC</a></span>(<span class='no'>f</span>[, <span class='st'>"FOCUS C"</span>])</div><div class='output co'>#> df BIC +#> SFO 3 59.88504 +#> FOMC 4 45.47542 +#> DFOP 5 30.00984</div><div class='input'> </div></pre> </div> 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 @@ <meta property="og:title" content="Confidence intervals for parameters of mkinfit objects — confint.mkinfit" /> -<meta property="og:description" content="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)." /> +<meta property="og:description" content="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)." /> <meta name="twitter:card" content="summary" /> @@ -71,7 +76,7 @@ could likely be improved by using the method of Venzon and Moolgavkar (1988)." / </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.6</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> @@ -135,14 +140,19 @@ could likely be improved by using the method of Venzon and Moolgavkar (1988)." / </div> <div class="ref-description"> - <p>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).</p> + <p>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).</p> </div> <pre class="usage"><span class='co'># S3 method for mkinfit</span> <span class='fu'><a href='https://rdrr.io/r/stats/confint.html'>confint</a></span>(<span class='no'>object</span>, <span class='no'>parm</span>, <span class='kw'>level</span> <span class='kw'>=</span> <span class='fl'>0.95</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>1</span> - - <span class='no'>level</span>, <span class='no'>cutoff</span>, <span class='kw'>method</span> <span class='kw'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span>(<span class='st'>"profile"</span>, <span class='st'>"quadratic"</span>), + <span class='no'>level</span>, <span class='no'>cutoff</span>, <span class='kw'>method</span> <span class='kw'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span>(<span class='st'>"quadratic"</span>, <span class='st'>"profile"</span>), <span class='kw'>transformed</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>backtransform</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>cores</span> <span class='kw'>=</span> <span class='fu'><a href='https://rdrr.io/r/base/Round.html'>round</a></span>(<span class='fu'>detectCores</span>()/<span class='fl'>2</span>), <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>, <span class='no'>...</span>)</pre> @@ -174,11 +184,11 @@ cutoff value overrides arguments 'level' and 'alpha'</p></td> </tr> <tr> <th>method</th> - <td><p>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.</p></td> + <td><p>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.</p></td> </tr> <tr> <th>transformed</th> @@ -231,20 +241,27 @@ machines, cores > 1 is not supported.</p></td> #> parent_0 73.0641834 92.1392181 #> k_parent_sink 0.2170293 0.4235348 #> sigma 3.1307772 8.0628314</div><div class='input'> +<span class='co'># Set the number of cores for the profiling method for further examples</span> +<span class='kw'>if</span> (<span class='fu'><a href='https://rdrr.io/r/base/identical.html'>identical</a></span>(<span class='fu'><a href='https://rdrr.io/r/base/Sys.getenv.html'>Sys.getenv</a></span>(<span class='st'>"NOT_CRAN"</span>), <span class='st'>"true"</span>)) { + <span class='no'>n_cores</span> <span class='kw'><-</span> <span class='kw pkg'>parallel</span><span class='kw ns'>::</span><span class='fu'><a href='https://rdrr.io/r/parallel/detectCores.html'>detectCores</a></span>() - <span class='fl'>1</span> +} <span class='kw'>else</span> { + <span class='no'>n_cores</span> <span class='kw'><-</span> <span class='fl'>1</span> +} +<span class='kw'>if</span> (<span class='fu'><a href='https://rdrr.io/r/base/Sys.getenv.html'>Sys.getenv</a></span>(<span class='st'>"TRAVIS"</span>) <span class='kw'>!=</span> <span class='st'>""</span>) <span class='no'>n_cores</span> <span class='kw'>=</span> <span class='fl'>1</span> +<span class='kw'>if</span> (<span class='fu'><a href='https://rdrr.io/r/base/Sys.info.html'>Sys.info</a></span>()[<span class='st'>"sysname"</span>] <span class='kw'>==</span> <span class='st'>"Windows"</span>) <span class='no'>n_cores</span> <span class='kw'>=</span> <span class='fl'>1</span> + <span class='no'>SFO_SFO</span> <span class='kw'><-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span>(<span class='kw'>parent</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>, <span class='st'>"m1"</span>), <span class='kw'>m1</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>), <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) <span class='no'>SFO_SFO.ff</span> <span class='kw'><-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span>(<span class='kw'>parent</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>, <span class='st'>"m1"</span>), <span class='kw'>m1</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>), <span class='kw'>use_of_ff</span> <span class='kw'>=</span> <span class='st'>"max"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) <span class='no'>f_d_1</span> <span class='kw'><-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='no'>SFO_SFO</span>, <span class='fu'><a href='https://rdrr.io/r/base/subset.html'>subset</a></span>(<span class='no'>FOCUS_2006_D</span>, <span class='no'>value</span> <span class='kw'>!=</span> <span class='fl'>0</span>), <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'><a href='https://rdrr.io/r/base/system.time.html'>system.time</a></span>(<span class='no'>ci_profile</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/stats/confint.html'>confint</a></span>(<span class='no'>f_d_1</span>, <span class='kw'>cores</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>))</div><div class='output co'>#> User System verstrichen -#> 51.646 0.000 51.673 </div><div class='input'><span class='co'># The following does not save much time, as parent_0 takes up most of the time</span> -<span class='co'># system.time(ci_profile <- confint(f_d_1, cores = 5))</span> -<span class='co'># system.time(ci_profile <- confint(f_d_1,</span> -<span class='co'># c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 1))</span> -<span class='co'># If we exclude parent_0 (the confidence of which is often of minor interest), we get a nice</span> -<span class='co'># performance improvement from about 30 seconds to about 12 seconds</span> -<span class='co'># system.time(ci_profile_no_parent_0 <- confint(f_d_1, </span> -<span class='co'># c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = 4))</span> -<span class='no'>ci_profile</span></div><div class='output co'>#> 2.5% 97.5% +<span class='fu'><a href='https://rdrr.io/r/base/system.time.html'>system.time</a></span>(<span class='no'>ci_profile</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/stats/confint.html'>confint</a></span>(<span class='no'>f_d_1</span>, <span class='kw'>method</span> <span class='kw'>=</span> <span class='st'>"profile"</span>, <span class='kw'>cores</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>))</div><div class='output co'>#> User System verstrichen +#> 51.297 0.000 51.328 </div><div class='input'><span class='co'># Using more cores does not save much time here, as parent_0 takes up most of the time</span> +<span class='co'># If we additionally exclude parent_0 (the confidence of which is often of</span> +<span class='co'># minor interest), we get a nice performance improvement from about 50</span> +<span class='co'># seconds to about 12 seconds if we use at least four cores</span> +<span class='fu'><a href='https://rdrr.io/r/base/system.time.html'>system.time</a></span>(<span class='no'>ci_profile_no_parent_0</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/stats/confint.html'>confint</a></span>(<span class='no'>f_d_1</span>, <span class='kw'>method</span> <span class='kw'>=</span> <span class='st'>"profile"</span>, + <span class='fu'><a href='https://rdrr.io/r/base/c.html'>c</a></span>(<span class='st'>"k_parent_sink"</span>, <span class='st'>"k_parent_m1"</span>, <span class='st'>"k_m1_sink"</span>, <span class='st'>"sigma"</span>), <span class='kw'>cores</span> <span class='kw'>=</span> <span class='no'>n_cores</span>))</div><div class='output co'>#> <span class='message'>Profiling the likelihood</span></div><div class='output co'>#> User System verstrichen +#> 0.006 0.003 11.435 </div><div class='input'><span class='no'>ci_profile</span></div><div class='output co'>#> 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.</p></td> <span class='co'># interval based on the untransformed fit for k_m1_sink</span> <span class='no'>rel_diffs_transformed</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/base/MathFun.html'>abs</a></span>((<span class='no'>ci_quadratic_transformed</span> - <span class='no'>ci_profile</span>)/<span class='no'>ci_profile</span>) <span class='no'>rel_diffs_untransformed</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/base/MathFun.html'>abs</a></span>((<span class='no'>ci_quadratic_untransformed</span> - <span class='no'>ci_profile</span>)/<span class='no'>ci_profile</span>) -<span class='no'>rel_diffs_transformed</span></div><div class='output co'>#> 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</div><div class='input'><span class='no'>rel_diffs_untransformed</span></div><div class='output co'>#> 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</div><div class='input'> -<span class='co'># Set the number of cores for further examples</span> -<span class='kw'>if</span> (<span class='fu'><a href='https://rdrr.io/r/base/identical.html'>identical</a></span>(<span class='fu'><a href='https://rdrr.io/r/base/Sys.getenv.html'>Sys.getenv</a></span>(<span class='st'>"NOT_CRAN"</span>), <span class='st'>"true"</span>)) { - <span class='no'>n_cores</span> <span class='kw'><-</span> <span class='kw pkg'>parallel</span><span class='kw ns'>::</span><span class='fu'><a href='https://rdrr.io/r/parallel/detectCores.html'>detectCores</a></span>() - <span class='fl'>1</span> -} <span class='kw'>else</span> { - <span class='no'>n_cores</span> <span class='kw'><-</span> <span class='fl'>1</span> -} -<span class='kw'>if</span> (<span class='fu'><a href='https://rdrr.io/r/base/Sys.getenv.html'>Sys.getenv</a></span>(<span class='st'>"TRAVIS"</span>) <span class='kw'>!=</span> <span class='st'>""</span>) <span class='no'>n_cores</span> <span class='kw'>=</span> <span class='fl'>1</span> -<span class='kw'>if</span> (<span class='fu'><a href='https://rdrr.io/r/base/Sys.info.html'>Sys.info</a></span>()[<span class='st'>"sysname"</span>] <span class='kw'>==</span> <span class='st'>"Windows"</span>) <span class='no'>n_cores</span> <span class='kw'>=</span> <span class='fl'>1</span> +<span class='no'>rel_diffs_transformed</span> <span class='kw'><</span> <span class='no'>rel_diffs_untransformed</span></div><div class='output co'>#> 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</div><div class='input'><span class='fu'><a href='https://rdrr.io/r/base/Round.html'>signif</a></span>(<span class='no'>rel_diffs_transformed</span>, <span class='fl'>3</span>)</div><div class='output co'>#> 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</div><div class='input'><span class='fu'><a href='https://rdrr.io/r/base/Round.html'>signif</a></span>(<span class='no'>rel_diffs_untransformed</span>, <span class='fl'>3</span>)</div><div class='output co'>#> 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</div><div class='input'> <span class='co'># Investigate a case with formation fractions</span> <span class='no'>f_d_2</span> <span class='kw'><-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='no'>SFO_SFO.ff</span>, <span class='fu'><a href='https://rdrr.io/r/base/subset.html'>subset</a></span>(<span class='no'>FOCUS_2006_D</span>, <span class='no'>value</span> <span class='kw'>!=</span> <span class='fl'>0</span>), <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='no'>ci_profile_ff</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/stats/confint.html'>confint</a></span>(<span class='no'>f_d_2</span>, <span class='kw'>cores</span> <span class='kw'>=</span> <span class='no'>n_cores</span>)</div><div class='output co'>#> <span class='message'>Profiling the likelihood</span></div><div class='input'><span class='no'>ci_profile_ff</span></div><div class='output co'>#> 2.5% 97.5% +<span class='no'>ci_profile_ff</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/stats/confint.html'>confint</a></span>(<span class='no'>f_d_2</span>, <span class='kw'>method</span> <span class='kw'>=</span> <span class='st'>"profile"</span>, <span class='kw'>cores</span> <span class='kw'>=</span> <span class='no'>n_cores</span>)</div><div class='output co'>#> <span class='message'>Profiling the likelihood</span></div><div class='input'><span class='no'>ci_profile_ff</span></div><div class='output co'>#> 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.</p></td> <span class='no'>rel_diffs_untransformed_ff</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/base/MathFun.html'>abs</a></span>((<span class='no'>ci_quadratic_untransformed_ff</span> - <span class='no'>ci_profile_ff</span>)/<span class='no'>ci_profile_ff</span>) <span class='co'># While the confidence interval for the parent rate constant is closer to</span> <span class='co'># the profile based interval when using the internal parameter</span> -<span class='co'># transformation, the intervals for the other parameters are 'better</span> +<span class='co'># transformation, the interval for the metabolite rate constant is 'better</span> <span class='co'># without internal parameter transformation.</span> -<span class='no'>rel_diffs_transformed_ff</span></div><div class='output co'>#> 2.5% 97.5% +<span class='no'>rel_diffs_transformed_ff</span> <span class='kw'><</span> <span class='no'>rel_diffs_untransformed_ff</span></div><div class='output co'>#> 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</div><div class='input'><span class='no'>rel_diffs_transformed_ff</span></div><div class='output co'>#> 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 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> 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." /> </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.6</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> @@ -167,12 +167,29 @@ as used internally during the optimisation?</p></td> <p>A numeric vector of fitted model parameters</p> + <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> + <pre class="examples"><div class='input'><span class='no'>fit</span> <span class='kw'><-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='st'>"SFO"</span>, <span class='no'>FOCUS_2006_C</span>)</div><div class='output co'>#> <span class='message'>Ordinary least squares optimisation</span></div><div class='output co'>#> 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</div><div class='output co'>#> <span class='message'>Optimisation successfully terminated.</span></div><div class='input'><span class='fu'>parms</span>(<span class='no'>fit</span>)</div><div class='output co'>#> parent_0 k_parent_sink sigma +#> 82.4921598 0.3060633 4.6730124 </div><div class='input'><span class='fu'>parms</span>(<span class='no'>fit</span>, <span class='kw'>transformed</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)</div><div class='output co'>#> parent_0 log_k_parent_sink sigma +#> 82.492160 -1.183963 4.673012 </div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> <h2>Contents</h2> <ul class="nav nav-pills nav-stacked"> <li><a href="#arguments">Arguments</a></li> <li><a href="#value">Value</a></li> + <li><a href="#examples">Examples</a></li> </ul> </div> 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 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.6</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.7</span> </span> </div> @@ -144,7 +144,7 @@ <colgroup><col class="name" /><col class="desc" /></colgroup> <tr> <th>object</th> - <td><p>An <code><a href='mkinfit.html'>mkinfit</a></code> object</p></td> + <td><p>A <code><a href='mkinfit.html'>mkinfit</a></code> object</p></td> </tr> <tr> <th>standardized</th> diff --git a/docs/reference/update.mkinfit-1.png b/docs/reference/update.mkinfit-1.png Binary files differnew file mode 100644 index 00000000..4cbefa56 --- /dev/null +++ b/docs/reference/update.mkinfit-1.png diff --git a/docs/reference/update.mkinfit-2.png b/docs/reference/update.mkinfit-2.png Binary files differnew file mode 100644 index 00000000..f432f6f8 --- /dev/null +++ b/docs/reference/update.mkinfit-2.png 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." /> </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">mkin</a> - <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.6</span> + <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.9.49.8</span> </span> </div> @@ -167,349 +167,11 @@ remove arguments given in the original call</p></td> <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> <pre class="examples"><div class='input'><span class='co'># \dontrun{</span> -<span class='no'>fit</span> <span class='kw'><-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='st'>"DFOP"</span>, <span class='fu'><a href='https://rdrr.io/r/base/subset.html'>subset</a></span>(<span class='no'>FOCUS_2006_D</span>, <span class='no'>value</span> <span class='kw'>!=</span> <span class='fl'>0</span>), <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'><a href='https://rdrr.io/r/stats/update.html'>update</a></span>(<span class='no'>fit</span>, <span class='kw'>error_model</span> <span class='kw'>=</span> <span class='st'>"tc"</span>)</div><div class='output co'>#> $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" </div><div class='input'># } +<span class='no'>fit</span> <span class='kw'><-</span> <span class='fu'><a href='mkinfit.html'>mkinfit</a></span>(<span class='st'>"SFO"</span>, <span class='fu'><a href='https://rdrr.io/r/base/subset.html'>subset</a></span>(<span class='no'>FOCUS_2006_D</span>, <span class='no'>value</span> <span class='kw'>!=</span> <span class='fl'>0</span>), <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) +<span class='fu'><a href='parms.html'>parms</a></span>(<span class='no'>fit</span>)</div><div class='output co'>#> parent_0 k_parent_sink sigma +#> 99.44423886 0.09793574 3.39632469 </div><div class='input'><span class='fu'><a href='plot.mkinfit.html'>plot_err</a></span>(<span class='no'>fit</span>)</div><div class='img'><img src='update.mkinfit-1.png' alt='' width='700' height='433' /></div><div class='input'><span class='no'>fit_2</span> <span class='kw'><-</span> <span class='fu'><a href='https://rdrr.io/r/stats/update.html'>update</a></span>(<span class='no'>fit</span>, <span class='kw'>error_model</span> <span class='kw'>=</span> <span class='st'>"tc"</span>) +<span class='fu'><a href='parms.html'>parms</a></span>(<span class='no'>fit_2</span>)</div><div class='output co'>#> parent_0 k_parent_sink sigma_low rsd_high +#> 1.008549e+02 1.005665e-01 3.752222e-03 6.763434e-02 </div><div class='input'><span class='fu'><a href='plot.mkinfit.html'>plot_err</a></span>(<span class='no'>fit_2</span>)</div><div class='img'><img src='update.mkinfit-2.png' alt='' width='700' height='433' /></div><div class='input'># } </div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> diff --git a/man/AIC.mmkin.Rd b/man/AIC.mmkin.Rd index a49b69b8..a10d0aeb 100644 --- a/man/AIC.mmkin.Rd +++ b/man/AIC.mmkin.Rd @@ -31,15 +31,21 @@ same dataset. 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"]) } } diff --git a/man/confint.mkinfit.Rd b/man/confint.mkinfit.Rd index ee07c9c1..e4a60556 100644 --- a/man/confint.mkinfit.Rd +++ b/man/confint.mkinfit.Rd @@ -5,7 +5,7 @@ \title{Confidence intervals for parameters of mkinfit objects} \usage{ \method{confint}{mkinfit}(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, ...) } @@ -23,11 +23,11 @@ confidence intervals. If missing, all parameters are considered.} in the log-likelihoods at the confidence boundary. Specifying an explicit cutoff value overrides arguments 'level' and 'alpha'} -\item{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.} +\item{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.} \item{transformed}{If the quadratic approximation is used, should it be applied to the likelihood based on the transformed parameters?} @@ -49,9 +49,14 @@ A matrix with columns giving lower and upper confidence limits for each parameter. } \description{ -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). } \examples{ f <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE) @@ -60,19 +65,26 @@ confint(f, method = "quadratic") \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 @@ -84,21 +96,14 @@ ci_quadratic_untransformed # 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 @@ -108,8 +113,9 @@ rel_diffs_transformed_ff <- abs((ci_quadratic_transformed_ff - ci_profile_ff)/ci 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 diff --git a/man/parms.Rd b/man/parms.Rd index 73cb23cd..5752de0b 100644 --- a/man/parms.Rd +++ b/man/parms.Rd @@ -25,3 +25,8 @@ 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. } +\examples{ +fit <- mkinfit("SFO", FOCUS_2006_C) +parms(fit) +parms(fit, transformed = TRUE) +} diff --git a/man/residuals.mkinfit.Rd b/man/residuals.mkinfit.Rd index 407b89b9..aaff12c0 100644 --- a/man/residuals.mkinfit.Rd +++ b/man/residuals.mkinfit.Rd @@ -7,7 +7,7 @@ \method{residuals}{mkinfit}(object, standardized = FALSE, ...) } \arguments{ -\item{object}{An \code{\link{mkinfit}} object} +\item{object}{A \code{\link{mkinfit}} object} \item{standardized}{Should the residuals be standardized by dividing by the standard deviation obtained from the fitted error model?} diff --git a/man/update.mkinfit.Rd b/man/update.mkinfit.Rd index aae1fbb4..7054d2e6 100644 --- a/man/update.mkinfit.Rd +++ b/man/update.mkinfit.Rd @@ -23,7 +23,11 @@ override these starting values. } \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) } } |
