Make the 'quadratic' the default for 'confint'

Also the documentation was improved here and there

5 files changed, 54 insertions, 34 deletions

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)