aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2019-10-28 16:39:14 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2019-10-28 18:32:00 +0100
commit900790b4139dd672c7383a3ed6ad2c1e51d855b9 (patch)
treecb77959bd7ce66b33c81d28676fc5ee87ae0bc4a /R
parentcc53cf26628a0433e6edd157c87edab340cdd013 (diff)
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
Diffstat (limited to 'R')
-rw-r--r--R/confint.mkinfit.R100
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), "%")

Contact - Imprint