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 --- docs/reference/confint.mkinfit.html | 109 +++++++++++++++++++++--------------- 1 file changed, 64 insertions(+), 45 deletions(-) (limited to 'docs/reference/confint.mkinfit.html') 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 -- cgit v1.2.1