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 --- man/AIC.mmkin.Rd | 8 +++++- man/confint.mkinfit.Rd | 66 ++++++++++++++++++++++++++---------------------- man/parms.Rd | 5 ++++ man/residuals.mkinfit.Rd | 2 +- man/update.mkinfit.Rd | 8 ++++-- 5 files changed, 55 insertions(+), 34 deletions(-) (limited to 'man') 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) } } -- cgit v1.2.1