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 --- R/confint.mkinfit.R | 64 +++++++++++++++++++++++++++++------------------------ 1 file changed, 35 insertions(+), 29 deletions(-) (limited to 'R/confint.mkinfit.R') 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, ...) { -- cgit v1.2.1