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/AIC.mmkin.R | 8 ++++++- R/confint.mkinfit.R | 64 ++++++++++++++++++++++++++++----------------------- R/parms.mkinfit.R | 6 ++++- R/residuals.mkinfit.R | 2 +- R/update.mkinfit.R | 8 +++++-- 5 files changed, 54 insertions(+), 34 deletions(-) (limited to 'R') 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) -- cgit v1.2.1