Parallel computation for confidence intervals

Only on Linux at the moment. Some more examples in the help page.
Remove the distribution argument for the quadratic method

1 files changed, 82 insertions, 18 deletions

diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R
index 8467a85b..75813360 100644
--- a/R/confint.mkinfit.R
+++ b/R/confint.mkinfit.R
@@ -22,15 +22,18 @@
 #' @param backtransform If we approximate the likelihood in terms of the
 #'   transformed parameters, should we backtransform the parameters with
 #'   their confidence intervals?
-#' @param distribution For the quadratic approximation, should we use
-#'   the student t distribution or assume normal distribution for
-#'   the parameter estimate
-#' @param quiet Should we suppress messages?
+#' @param cores The number of cores to be used for multicore processing. This
+#'   is only used when the \code{cluster} argument is \code{NULL}. On Windows
+#'   machines, cores > 1 is not supported.
+#' @param quiet Should we suppress the message "Profiling the likelihood"
 #' @return A matrix with columns giving lower and upper confidence limits for
 #'   each parameter.
 #' @param \dots Not used
 #' @importFrom stats qnorm
-#' @references Pawitan Y (2013) In all likelihood - Statistical modelling and
+#' @references
+#'   Bates DM and Watts GW (1988) Nonlinear regression analysis & its applications
+#'
+#'   Pawitan Y (2013) In all likelihood - Statistical modelling and
 #'   inference using likelihood. Clarendon Press, Oxford.
 #'
 #'   Venzon DJ and Moolgavkar SH (1988) A Method for Computing
@@ -39,15 +42,78 @@
 #' @examples
 #' f <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE)
 #' confint(f, method = "quadratic")
+#'
 #' \dontrun{
-#'   confint(f, method = "profile")
+#' confint(f, method = "profile")
+#'
+#' 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))
+#' ci_profile
+#' ci_quadratic_transformed <- confint(f_d_1, method = "quadratic")
+#' ci_quadratic_transformed
+#' ci_quadratic_untransformed <- confint(f_d_1, method = "quadratic", transformed = FALSE)
+#' ci_quadratic_untransformed
+#' # Against the expectation based on Bates and Watts (1988), the confidence
+#' # intervals based on the internal parameter transformation are less
+#' # congruent with the likelihood based intervals. Note the superiority of the
+#' # 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
+#'
+#' # 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
+#' ci_quadratic_transformed_ff <- confint(f_d_2, method = "quadratic")
+#' ci_quadratic_transformed_ff
+#' ci_quadratic_untransformed_ff <- confint(f_d_2, method = "quadratic", transformed = FALSE)
+#' ci_quadratic_untransformed_ff
+#' rel_diffs_transformed_ff <- abs((ci_quadratic_transformed_ff - ci_profile_ff)/ci_profile_ff)
+#' 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
+#' # without internal parameter transformation.
+#' rel_diffs_transformed_ff
+#' rel_diffs_untransformed_ff
+#'
+#' # The profiling for the following fit does not finish in a reasonable time
+#' #m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")),
+#' #  M1 = mkinsub("SFO"),
+#' #  M2 = mkinsub("SFO"),
+#' #  use_of_ff = "max", quiet = TRUE)
+#' #DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data
+#' #f_tc_2 <- mkinfit(m_synth_DFOP_par, DFOP_par_c, error_model = "tc",
+#' #  error_model_algorithm = "direct", quiet = TRUE)
+#' #confint(f_tc_2, "parent_0")
 #' }
 #' @export
 confint.mkinfit <- function(object, parm,
   level = 0.95, alpha = 1 - level, cutoff,
   method = c("profile", "quadratic"),
   transformed = TRUE, backtransform = TRUE,
-  distribution = c("student_t", "normal"), quiet = FALSE, ...)
+  cores = round(detectCores()/2), quiet = FALSE, ...)
 {
   tparms <- parms(object, transformed = TRUE)
   bparms <- parms(object, transformed = FALSE)
@@ -68,11 +134,7 @@ confint.mkinfit <- function(object, parm,
 
   if (method == "quadratic") {
 
-    distribution <- match.arg(distribution)
-
-    quantiles <- switch(distribution,
-      student_t = qt(a, object$df.residual),
-      normal = qnorm(a))
+    quantiles <- qt(a, object$df.residual)
 
     covar_pnames <- if (missing(parm)) {
       if (transformed) tpnames else bpnames
@@ -99,7 +161,7 @@ confint.mkinfit <- function(object, parm,
       ses     <- sqrt(diag(covar))[covar_pnames]
       lci    <- covar_parms + quantiles[1] * ses
       uci    <- covar_parms + quantiles[2] * ses
-      if (backtransform) {
+      if (transformed & backtransform) {
         lci_back <- backtransform_odeparms(lci,
           object$mkinmod, object$transform_rates, object$transform_fractions)
         lci <- c(lci_back, lci[names(object$errparms)])
@@ -108,6 +170,7 @@ confint.mkinfit <- function(object, parm,
         uci <- c(uci_back, uci[names(object$errparms)])
       }
     }
+    ci <- cbind(lower = lci, upper = uci)
   }
 
   if (method == "profile") {
@@ -125,8 +188,7 @@ confint.mkinfit <- function(object, parm,
 
     all_parms <- parms(object)
 
-    for (pname in profile_pnames)
-    {
+    get_ci <- function(pname) {
       pnames_free <- setdiff(names(all_parms), pname)
       profile_ll <- function(x)
       {
@@ -143,12 +205,14 @@ confint.mkinfit <- function(object, parm,
         (cutoff - (object$logLik - profile_ll(x)))^2
       }
 
-      lci[pname] <- optimize(cost, lower = 0, upper = all_parms[pname])$minimum
-      uci[pname] <- optimize(cost, lower = all_parms[pname], upper = 15 * all_parms[pname])$minimum
+      lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname])$minimum
+      uci_pname <- optimize(cost, lower = all_parms[pname],
+        upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]))$minimum
+      return(c(lci_pname, uci_pname))
     }
+    ci <- t(parallel::mcmapply(get_ci, profile_pnames, mc.cores = cores))
   }
 
-  ci <- cbind(lower = lci, upper = uci)
   colnames(ci) <- paste0(
     format(100 * a, trim = TRUE, scientific = FALSE, digits = 3), "%")